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THE DYNAMO MAINTENANCE OF STEADY 
MAGNETIC FIELDS 


By T. G. COWLING (The Unive rsity, Leeds) 


teceived 5 July 1956] 


SUMMARY 


A general type of mathematical argument is described, which applies to all the 
eases in which dynamo maintenance of a steady magnetic field by motion in a 
uniform mass is known to be impossible. The number of general classes of motion 
to which this type of argument applies is shown to be strictly limited. It is suggested 
that all proofs of the impossibility of dynamo maintenance for general classes of 
motion must be reducible to this type. 


1. THE steady-state dynamo problem in a uniform mass can be stated as 
follows. Given a uniform mass of electrically conducting liquid, contained 
ina volume 7), what condition must be satisfied by the velocity v of a 
steady motion in 7, in order that a steady magnetic field H can be main- 
tained by the dynamo interaction between the motion and the field? By 
numerical solution of the differential equations (1), maintenance has been 
shown always to be possible for certain flow patterns if the flow is sufficiently 
fast; that is, the maintenance can always be secured by increasing a given 
initial flow velocity v in a suitable ratio A:1. For certain other flow patterns 
dynamo maintenance has been shown to be impossible, however much the 
velocity Vv is increased (2, 3,4). Such flow patterns include two-dimensional 
motions, motions with axial symmetry, and ‘toroidal’ motions (in which 
the stream-lines lie on concentric spheres). 


A closer study of the cases in which dynamo maintenance of a steady 


field is impossible may throw light on the general dynamo problem. Such 


a study is attempted here. The results suggest that the number of cases 
in which dynamo maintenance is impossible is likely to be a very restricted 


group. 


2. A transformation of the dynamo equations is first given for a non- 
steady field (cf. (5)). If E, H are the electric and magnetic intensities,+ 
the electric force on the (non-magnetic) material is E+-v , H; thus if j is 
the current density, o the electrical conductivity, 


j — o(E+vH). (1) 


Electromagnetic units are used throughout. 
[Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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, 
Let H’ be a second field, corresponding to an arbitrary current system j’ 
within the volume 7, of the conducting liquid. Then 
eo r oie i 
j.ji' dz, (E+vAH).j’ dry. (2) 
Co ‘ 
The integrals in this equation can, if desired, be supposed to be taken over 
the whole volume of space, since j’ = 0 outside 7); when this is done, the 
volume element dz, will be replaced by dr. Since 
curl H’ = 47j’ (3) 
it follows that 
, 2: 9 lf. ae , , >) 1 
E.j' dz, | E.curl H’ dr fdiv(H’ , E)+H’.curl E} dr. 
, di p, 4a an 
(4) 


Now, by Green's theorem, 
div(H’ A E) dr 0. 


The surface integral over an infinite sphere introduces no complications, ' 
since H’ and E vanish at large distances r like r~-* and r~?; the boundary 
of 7, equally introduces no difficulties, since the normal component of 
H’ \ E depends only on the tangential components of H’ and E, and so 


is continuous at the boundary. Again. 


. H 
curl E : (5) 
ot ' 
and so equation (4) becomes 
r - 1 (,,, 0H 
| E.j’ dr | w= de. 
F 4a. ot 
Using this in equation (2) we have finally 
ree ‘ - : fas, OE 
| j.j' dz, | (v\H).j’ dz, - | n°. = dr. (6) 
o. J dr | ct 


This equation, with arbitrary j’, is equivalent to the differential equation 

for H, and provides the desired transformation of the dynamo equations. 
3. Now suppose that j’ can be found such that 

(vA H).j' = v.(H/j’) = v.grady, (7) 

where us is some scalar function. Then, since v satisfies an equation of 

continuity dine 0. 


it follows that [ (vA H).j’ dz, ( div(ves) dro. 





T 
Thi: 


vani 


Ace 


The 
if th 
left « 


used 


4, © 

Su 
and ‘ 
Viss 
(6) h 
toen 
syste 
sectic 

Th 
perpe 
from 
equat 


Hene 


readi 


This 

H. 
Su} 

from - 


Take 
Satisfi 
vanis} 














THE DYNAMO MAINTENANCE OF STEADY MAGNETIC FIELDS 131 
| This vanishes, by Green’s theorem, since the normal component of v must 


Thus, in this case 


[ ii’ ar, eet 
G 3 ct 


« 


vanish on the boundary ol Tt 


0° 


Ir. 
Accordingly, the motion can maintain a steady field only if 


| ii dry = 0. (8) 


The maintenance of a steady field can therefore be proved to be impossible 
if the j’ satisfying equation (7) can be shown to make the integral on the 
left of equation (8) essentially positive. An argument of this type can be 
used for each of the cases mentioned in section 1. 


4, Case I: two-dimensional motion 
Suppose that 7, is an infinite cylinder with generators parallel to Oz, 

ind that v and H have cartesian components depending only on 2 and y; 

v is supposed to be perpendicular to Oz. The volume integrals in equation 

f (6) have then to be limited to the region between two planes z = const; 

} toensure the convergence of the last of the integrals at infinity, the current 


system j’ is assumed to be such as to give zero total current across any 


SO section of the cylinder. 
The field H can be divided into two parts, respectively parallel and 
perpendicular to Oz. The first, which vanishes outside the cylinder, arises 
| from the part of j perpendicular to Oz. Take this part of j as j’; then an 
equation of form (7) is found to be satisfied, with 


ob i O7r. 


Hence the first term on the right of equation (6) vanishes; the equation is 
readily found to reduce to 


Pf a l . ‘ 
H? dr, grad H.|* dtp. (9) 
6 ct . 270 | ' 
This shows that no steady field is possible unless grad H, = 0, giving 
tion H. = 0. 
ions. Suppose now that H, = 0. The field perpendicular to Oz can be derived 
irom the z-component of the vector potential A by the formulae 
H, cA,/cy, H, CA,/ 0x. 


Take j’ equal to this z-component of A; then an equation of type (7) is 
on of satisfied, with u& $A>. Hence the first integral on the right of (6 
vanishes; also, using Green’s theorem 


| . . 


— ] 
j.j' dr, V?A..A, dr = grad A, |? dr. 
4a — 4 . 


- . 7 


y 
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“ee, by a similar transformation, 





, CH aoe ¥ , OA, 
HY, — dr: | H’.—curl A, dr = | curl H’.? dr = 2a — o [4 dry. 
ct J ct P : t at | 
Thus equation (6) becomes 
or wz Es, ‘ 
= | A? dry = —5— | grad A. |? dz. (10) 
This shows that maintenance of a steady field is impossible if A. + 0, so 
: ] z 


that H must be identically zero. 

It may be noted that this proof applies unchanged, save for an inversion 
of the order of the two steps, if v has a component parallel to Oz which, 
like the other components, depends only on x and y. 


5. Case II: axial symmetry 

Suppose that 7, and the vector fields v, H have axial symmetry about 02; 
the lines of the flow v are assumed to be in meridian planes through Oz. 
Take cylindrical polar coordinates R, 4, z. The fact that the currents j are 
confined to the interior of 7) implies that Hy, the ¢-component of H, 
vanishes outside and on the boundary of 75. 

Suppose that the field H’ is given by H,/R*. This field is singular on Oz, 
the current system j’ corresponding to it including a current J along Oz. In 
consequence, the argument leading to equation (6) needs modification. It 
can be shown that the resulting modification of equation (6) is simply that 
the left-hand side has to include an extra term 

] 


oO 


« 


Ij. dz 


the integral being taken along Oz. 
With this choice of H’, an equation of form (7) is satisfied, with 
us H* 87 R?. 


Also, the component j, of j makes no contribution to j.j’, so that 


Jd 





l 67° 


— arama 


¢ 
167? . ss R 





= e. ~ | ors an Aen am! Hg | dro. 
ren? | \\8"¢ | +e are) 


Since dr, = Rd Rdddz, on integration with respect to R and 4, 


dl | 2H c Hg dr, on | = de. 
1677? . R2 éR\R 827 J \ R] pao 
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Also on Oz, I (H,/2.R)p-. = 7,. Hence 


 . ee eo LF H,)\ |* l f¢ /A,\? 
| j.j' dr, | I}, dz _— | | gradi $\| dr, +— mg dz. 
d ‘ LG" J | \R}| 87 J \ FR} R-0 
Using this result, the modified form of equation (6) can be reduced to 
0 ( ; £ H,\? 1 ¢ /A,\? 
£ dry grad —¢| dr, —- | $ dz. (11) 
al 2r0 . R | a) \R/} p-o 
Thus no maintenance of a steady field is possible if Hy + 0. 

Suppose therefore that H, = 0. Then the vector potential A from which 
H is derived is wholly in the ¢-direction. Take j’ R°A; then an equation 
of the form (7) is satisfied, with » = }R?A*. Also, transforming by Green’s 
theorem, 

bP aa ee LP oe A 
| H’.—d- | H’. curl dr 
£77 ot ‘77 ct 
® A lo [ pogo 
| curl H’. — dr “ | R242 dr. 
T ct 2c. 
\gain 
, ia, 
j.j dz | VA. RPA dr 


L | |raveyRd)—24 (RA)\ dr. 
ir J | oR | 


The second term in the brackets { | goes out on integrating with respect 
R: thus 2 . 


Bs l " 
| j.j' dz, | orad( RA) |? dr. 
e a7 | 
Hence equation (6) now gives 

a ae l 7 > 9 5) 
| RA? dry | grad( RA)? dr, (12) 

ct ~ no 
owing that steady dynamo maintenance is impossible with A + 0. Thus 


steady field can exist. 

The argument applies without change, save for an inversion of the order 
t the two steps, if v has in addition a component in the ¢-direction, 
symmetric about Oz. 


6. Case III: toroidal motion 
Let the boundary of 7, be a sphere, centre O, and let 
v r A crad U, 


where r is the displacement from O, and U is an arbitrary function of 
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position. Then the motion is toroidal, the lines of flow being the inter- 
sections of surface U’ = const with spheres 7 = const. 
Let j’ be defined by 


4zj’ = curl(r(r.H)) —r/ grad(r.H) 


inside 7,, and j’ = 0 outside, so that j’ is toroidal. The corresponding field 
is given by H’ = r(r.H)+ grad Q inside 7), H’ = grad 2’ outside, where 
Q = Q’ on the boundary of 7, in order to make the tangential component 
of H’ continuous. A relation of form (7) is now satisfied with ys = (r.H)?/8z, 
so that the first term on the right of equation (6) vanishes as before. 

In the second term on the right of equation (6), the contribution from 
the parts grad Q, grad Q’ of H’ can readily be shown to vanish; thus the 
term becomes 
ve .00). 2 ar, i (fc) ds... 


ct Sar c 


. . 


An 
Also, using Green’s theorem and the relation div H 0, 
167? | j.j' dry | curlH.curl(r(r.H)) dz 


{ (r.H)r.curl curl H dz | (r. Hr. V2H dz 


( (r.H)V2(r.H) dr = | \grad(r.H))? dr. 
Hence equation (6) now becomes 


(r.H)? dz, grad(r.H)\? dz, 


27a | 


. 


c 


so that no steady field with a non-zero radial component is possible. 
We are therefore restricted to toroidal fields, of form 


H r A grad P. 


The field given by this equation is unaltered if an arbitrary function of r 
is added to P; thus P can be chosen so that 


( [ Pas —6 


integrated over the surface of all spheres centre O. Since H 0 outside 7p, 
P = 0 outside and on the boundary of 7,. Let H’ be a further toroidal 
field, such that j’ has a radial component P/r within 7,; this condition is 
sufficient to fix H’ uniquely. Then an equation of form (7) is satisfied, 
with w& }P?; thus the first integral on the right of equation (6) again 
vanishes. 


Again 
S 


eurl[r(r.j’)] = curl(r P) H. 
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er Hence, on again using Green’s theorem 
t H ee ij 
H’ .<— dr | H’.curl fe2 dr 
rr ct ior et 
a7 F es; ° eP 
curl H e.5 dr | P— dry. 
7 ct ct 
eld . } ° 
Also 
ere . | . 
ont j.j dr, = os j’.curlecurl(rP) dz, 
Q . a 
B. woes ig ; 
j’.{V?(rP)—grad div(rP)} dr, 
to 
om « 
the 


a ¥ 
= j’.{rV*(P)+grad|2P—div(rP)]! dzp. 
But since div j’ 0 
( j'.grad[2P—div(rP)] dz, { div{j’[2P—div(rP)}} dz, 0 


by Green’s theorem. Hence 


[ j.i' dz, — ( PV?P dr, —— | grad P|? dtp. 


Using these results, equation (6) now gives 


I ; 
P? dr | grad P |? dry. 
ct =71O0 


« . 


Thus a steady non-zero field (P 0) is impossible. 


7. One may ask to what other cases a similar argument is applicable. 
With a wide family of flows v, a j’ can easily be found such that an equation 
of form (7) is satisfied. For example, suppose that the lines of flow lie on 
the family of surfaces ¢ = const, so that v.grad¢@ = 0. Then an equation 


of form (7) is valid with j grad y A grad d, pb ty*, where y = H.gradd. 


g 
a) | ° ey 
It can also be shown that with this value of j’, 
( 
H F , 
 . dr 2a x” dry. 

ct a. 
‘ The first parts of the proofs given above for cases I, II, and III do, in fact, 
rest on the existence of a function ¢; in ease I, ¢ is z, in case II it is the 
lda 
a polar angle d, and in case IIT it is 37°. 
is 2 
‘od The difficulty comes when one considers the remaining integral in 
CU : 


equation (6), 


j J dtp. 


To establish a contradiction, and so prove dynamo maintenance impossible, 
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this integral had to be expressed as a positive or negative definite form, 


If, as above, + 
; j grad y \ grad¢ = curl(y grad ¢) 
then, provided that certain surface integrals vanish, it can readily be 
shown that , ir 

j.j dr = —- x V*H. grad ¢ dzp. 


+77 
ny 


Since y = H.grad 4, this integral can be shown to be definitely positive 
or negative only if, for any field H such that divH = 0, V?H.grad¢ can 
be expressed in terms of H.grad ¢; i.e. the component of V?H parallel to 
gradd depends only on the component of H in the same direction. This 
was true with each of the functions ¢ appropriate to cases I, II, and ITI; 
but such a property can clearly exist only for a very limited class of such 
functions. Thus the argument given for these cases cannot be generalized 
to apply to more than a limited class of types of motion. 


8. The suggestion appears plausible that the only complete proofs of the 
impossibility of steady-state dynamo maintenance must be of the type 
considered above. Equation (6), with arbitrary j’, is completely equivalent 
to the differential equation for H. Thus one would expect that dynamo 
maintenance can be proved impossible only if, for arbitrary H, we can find 
j/ such that (6) is impossible with CH/ct = 0. Since v can be increased in 
an arbitrary ratio, it is plausible to assume that the impossibility arises 
only when j’ can be found such that the first term on the right of (6) vanishes; 
since the only general restriction imposed on v is the equation div v = 0, 
this suggests that an equation of form (7) must be satisfied; and so on. 
The argument is certainly not complete, but it is plausible. If it is indeed 
sound, further search for a general theorem on the impossibility of dynamo 
maintenance would appear pointless. 
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MAGNETO-HYDRODYNAMICS OF A FINITE 
ROTATING DISK 


By K. STEWARTSON 


(Department of Mathematics, The University, Bristol) 
Received 17 May 1956] 


SUMMARY 


The effect of a vertical magnetic field on the motion of a shallow dish of mercury, 
t of whose base is rotating while the rest is fixed, is considered. It is found that 
in the absence of the field and when the field is large there is a steady solution 
which the mercury is either uniformly rotating or at rest, there being a thin 
tion layer separating the two regions. A comparison with experiment suggests 
hat the flow in the absence of the field is unstable but that the field stabilizes the 
‘tion while it is still not large enough to make substantial modifications to the 
city field in the friction layer. 


— 


Introduction 

Ix arecent paper Lehnert (1) has described a series of beautiful experiments 
in magneto-hydrodynamics. A shallow layer of mercury was placed on 
. horizontal plane, an annulus of which could rotate about its centre in its 
own plane, while the rest of the plane remained at rest. At the same time 
i vertical magnetic field H was imposed. It was observed that when H = 0 
the motion of the surface of the mercury consisted of a general rotation and 
it was practically plane. On the other hand, when H was large (of the order 
of 8,000 gauss) only the fluid above the annulus rotated, almost as if solid, 
while the rest of the fluid remained substantially at rest. 

It is proposed, in this paper, to offer a partial theoretical explanation of 
these phenomena. Attention is confined to the flow due to a finite horizontal 
copper disk of radius a rotating in its own plane (z = 0) while the rest of 
the plane is at rest. Above this plane is a layer of mercury of depth d and 
above that air. A vertical magnetic field H is imposed. For large values 
of H, which are considered first, a solution is obtained in which the fluid 
has the same angular velocity as the point on z = 0 vertically below it, 
except when r, the distance from the axis of symmetry, is nearly equal to a. 
In this region a thin friction layer is formed in which the azimuthal compo- 
nents of the magnetic field and the velocity change rapidly. The free 
surface is plane if ry > a and parabolic if r < a with a discontinuity in 


[Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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gradient at r =a. At first sight it appears that this solution agrees with 
Lehnert’s photographs with H = 8,000 gauss, provided we assume that 
the discontinuity in the gradient of the free surface at 7 = a would cause 
such ‘vortex rows’ as were observed by Lehnert. However, on considering 
the orders of magnitude of the neglected terms it follows that the theory 
can only apply at much larger values of H than those used by Lehnert and 
that even at the largest field he used the magnetic body forces were smaller 
than the inertia terms in the momentum equations. 

Accordingly in section 3 we consider the same problem when the magnetic 
field is absent. If the radius a of the disk were infinite the mercury would 
rotate as if solid while if a were zero it would be at rest everywhere. It is 
assumed, therefore, that except for a narrow layer near r = a these con- 
ditions also occur in the present problem, i.e. the fluid rotates almost as 
if solid ifr < aandis almost at rest ifr > a. The flow in the neighbourhood 
of r = a is then found and the hypothesis is shown to be consistent so far 
as the primary flow is concerned. Indeed, in the conditions of the ex- 
periment a = 25 em., d = 0-6 em., Q = 1, and the kinematic viscosity 
v = 1-2 10-3em.*sec.~!, so that the thickness of the narrow layer would 
be only about 0-1 cm. and the form of the free surface almost the same as 
when H were large. However, it is known that in problems involving 
rotating disks the secondary flow is often of importance in the determina- 
tion of the primary flow (2). It may be, therefore, that the hypothesis is 
incorrect, although it is difficult to give a satisfactory alternative 
picture. 

Since the solution in section 3 involves a narrow vorticity layer not 
attached to a boundary along its length it will almost certainly be unstable; 
and this would explain why Lehnert’s photographs when H = 0 do not 
agree with it. On the other hand, when H = 8,000 gauss and the magnetic 
body forces are still not large enough to dominate the inertia terms the 
photographs agree with the theory. Further it is known from the theoretical 
papers of Chandrasekhar (3), Stuart (4), and others that magnetic fields 
do tend in many cases to stabilize the motion of the fluid. It is inferred, 
therefore, that at the highest fields used by Lehnert the flow may be deduced 
from section 3, the effect of the magnetic field being principally to ensure 
stability. 

Lehnert (1) concluded from his experiments that the effect of the 
magnetic field was destabilizing because the flow with H large was markedly 
different from that with H zero, although in both cases the free surfaces 
looked steady. If the conclusions of this paper are correct, then, in spite 
of the smooth appearance of the surface, the interior of the fluid must be 
turbulent in the absence of the magnetic field. 
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2. Strong magnetic field 

The disk is of radius a and is rotating in its own horizontal plane about 
itscentre O with angular velocity Q. The rest of the plane of the disk (z = 0) 
is fixed. Let us choose cylindrical polar coordinates, r, 6, z, with origin at O, 
and with Oz as vertical axis of rotation; from symmetry all dependent 
variables are independent of @. Let d be the depth of mercury above the 
plane z = 0, H the magnitude of the imposed vertical magnetic field, and 
v,H the velocity and magnetic fields in the mercury. Assuming steady 


motion the equations are (5) 


V*H 470 curl(v \H) (2.1) 
(v.V)v (/477p)(H A curl H)—(1/p)grad P+vV2v (2.2) 
und div Vv div H 0, 2.3) 


where p is the permeability, o the conductivity, p the density, p (= P—gz) 
the pressure, g gravity, and v the kinematic viscosity of the mercury. The 
electric field E is then given by 
curl H toro(cE- pV AH), (2.4) 
where c is the velocity of light. 
First consider the boundary conditions on v. At z = 0 they are simply 


v 0 if a and v (0, Qr,0) if r <a. (2.5) 
If during the motion the equation of the free surface is given by 
d f(r) 

we shall assume that the motion is such that f’ is small. Then, if (w, v, w) 

re the components of v, and since the free surface always consists of the 
same particles, we have 

D(z—d—fi(r))/Dt 0, le.w=0 atz=d-f. (2.6) 

Further, the stress across the free surface must be continuous, and since 


outside the fluid it is atmospheric and normal to the free surface 
p = const., év/i 0, and ceuj/ez =O atz=—d-f. (2.7) 


Except for the first condition, it is legitimate to neglect f and apply them 
d and then to use the condition P const. to determine the equa 
tion of the free surface. 

For the boundary conditions on H we shall assume that the material 
below the plane z 0} is a very good conductor so that in it we may take 
©. Further above the free surface of the fluid there is a poor conductor 

such as air in which o may be set equal to zero. 
\t the boundary z () the tangential components of E and the normal 


component of H must be continuous. These conditions are not all inde- 
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pendent and reduce to two only. The discontinuities in the other com- 
ponents of E and H are accounted for by surface distributions of charge 
and current. On the upper side of z = 0 


1 (eH, a 





c 4g = 
4a | er Cz 


and on the lower side EH, = 0 since o = ©. 
Hence one boundary condition is that 
OH,/ér = H,|éz atz 0. (2.8) 
Similarly by comparing the values of cH, on the two sides and noting that 
H, is continuous, the second condition is that 
CH,/ez = 0 atz=0. (2.9) 
On the free surface all four tangential components of E and H must be 
continuous. In a region of zero conductivity curl E = curlH = 0 so that 
the tangential components of E and H must satisfy the irrotational con- 
dition at the free surface which may be taken to be z = d; 
1.e. E, = Hp = 0 


at z = d. (2.10) 


If the imposed magnetic field is strong enough a solution may be found 
on neglecting u, w, H,, and H,—H. 


Writing Ht = (0, 2, #), V (0, v, 0) (2.11) 
the equations reduce to 
steak he 1 ch . PO 
cz" cr” vr c7 i C2 
(2.12) 
Fara) Pate} map d pH ch 
cz? sor? ror r2 4izpv cz 
2H h* rs [P | 
dar p or Sz 
se sia (2.13) 
l h? 
and 0 : (P+ = 
p cz ST 
with boundary conditions 
—=h=0 atz=d, 
Ce 
2 Q) ch =( ifz Or<a (2.14) 
Cz 
] fe 
and e= et Oo t2—0¢> 6 
C2 





These equations are exact, but (2.13) leads to a contradiction since 


47v?—h?/p is a function of both r and z. When the imposed field is strong 








M 
enough 
a negli 
conside 
in orde 


(2.12) s 


where 
xy 


Furthe 


If His 
contril 
those f 
that th 
of that 
contro! 
than t] 
can be 
simplif 


°9 
how co- 


with b 


and 


The 


i 


where 














MAGNETO-HYDRODYNAMICS OF A 





FINITE ROTATING DISK 141 


enough, however, the values of w and w necessary to cancel this term have 
a negligible effect on (2.12). We proceed to solve (2.12) and then from a 
consideration of (2.13) we may find the condition which H must satisfy 
in order that the assumption may be justified. The general solution of 


2.12) satisfying the boundary conditions at z = d is 


y ( J,(kr)| A,(k) cosh «,(d—z)+ A,(k) cosh a,(d -z)| dk 
. , (2.15) 
h = 4n(pvo)? ( J,(kr)| A,(k) sinh a,(d—z)+ A,(k) sinh «,(d—z)] dk 


where A,, A, are unknown functions of k, 


B+ (kh §) Xo 3+-(k?+-B?)?, and 28 = wH(o/pv)*. (2.16) 
Further, since ch/éz 0 at 2 0, 
x, A, cosh x, d vy A, cosh vy d. (2.97) 


> 9 


If H is large so is f and therefore a, ~ 28, a, ~ k?/28. Thus from (2.17) the 
contributions to v and h from A, may be neglected in comparison with 
those from A, except for the contribution to éh/éz near z 0. This means 
that the solution in the main body of the fluid is practically independent 
fthat boundary condition. Moreover, since only the terms involving x, 
mtrol this solution the rate of variation in the z-direction is much smaller 


in that in the r-direction in the vicinity of r = a. Although a solution 


n be completed using (2.15) it is more convenient to use these results to 
plify the equations of motion. The dominating operators in (2.12) are 


w o*/er? and He/cz so that they reduce to 


2), ov on pl C h 


—4rucH and (2.18) 
( or tirpv C2 
th boundary conditions 
Qr. » v 0, r>a atz Vv 
0 atz=d 
The solution is 
LO ()) ( sin k(r—a cosh(k?/28)(d—z) di 
‘ - h cosh dk?/2B 
' ” ¥Iq\ 
: S { sn COSY (% : I} con k. (r -a)jexp|—k, |\r—a } +8, 
n=0 (2.19) 


where 
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while 
_—— sinf(n+4)x(1—z/d)} . 
h = 4Qr(pvo)! >| r (n # 21} sin k,|r—a\exp|[ —k,, |r—a)]. 
aie : (2.21) 
The boundary condition on h at z = 0 can now be satisfied by evaluating 
ch/éz at z = 0 from (2.21) and cancelling it with a boundary layer on z = 0. 
Thus if (2.21) gives a value 47(pvo)!BOQr g(r), say, for ch/éz at z h we add 
to (2.21) the terms 
h 27(pvo)!Qrg(r)e-2F> and v LOr g(r)e-2P2, (2.22) 


Since f is large the contribution from (2.22) to the boundary condition 
on vatz = Ois negligible so that (2.21) and (2.22) together give a uniformly 
valid solution of (2.12) and (2.14) as B > o. 

Of special interest is the velocity variation on the free surface where 
Qr Ss (—)" cos k,,(r- a)exp(—k, |r—a F 5 (2,23) 

Pa ; 2.2% 
0 


n--s 


7 


so that the change from a rigid body rotation to rest occurs in a distance 
O(d/B)'. In the strongest fields used by Lehnert (1) H = 8,000 gauss so that 
(d/8)* = 0-1em. Thus if this solution is appropriate outside a narrow layer 
of thickness a few millimetres the velocity is undisturbed, and through this 
layer the azimuthal velocity falls rapidly from a value Qa to zero. Although 
this is in substantial agreement with the photographs it will be shown 
below that the value 8,000 gauss for H is too low for the theory to apply 
and the magnetic forces are smaller than the others. However, the steady 
flow in the absence of a magnetic field is of the same type as this one with 
an even thinner friction layer and it is inferred that the magnetic field 
stabilizes the flow. 
The equation of the free surface is calculated from (2.13). Since 
4n(pvo)! = 5x 10-3 

we may neglect h? in comparison with v? and since the pressure is constant 
on it 


gz [ v?/r dr = const. 
0 
a ai ; i . 
rhus if a—r > (d/B)}, z= rt A, (2.24) 
2g 
where A is a constant nearly equal to d, and if r—a > (d/B)!, 
()2 ()2q2q} 
oe gt sso). (2.25) 
2g | gB* J 


In the intervening region the equation is more complicated, but the extra 
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terms are all small so that it is legitimate to neglect them and to regard the 
free surface as having a discontinuity in slope at r = a, being parabolic if 
a and horizontal if r a. 

The generalization to Lehnert’s experiments is now immediate, bearing 
in mind that the friction layer is confined to the vicinity of r = a, so that 
ifinstead of a disk there was a rotating annulus of inner and outer radii a 
ind b, the friction layers at 7 aandr b would be independent. Hence 
the fluid in between would rotate as if solid and the fluid outside them 
would be at rest. The free surface would consist of two planes joined by 
, paraboloid of revolution and the discontinuities in the tangents to the 
free surface at 7 = aandr = 6) would explain the distortion of the reflection 
f the grid in the free surface, which was observed by Lehnert. 

The principal assumptions we have made in the section are that the free 
surface is almost horizontal, i.e. |f’(r) 1 and that the inconsistency 
2.13) may be neglected. In the experiments 2 = 1 and a & 25 cm. so 


that at ost 2 5 
lat at m f Q*a/g =~ 0-025 


ind is small. The other difficulty is more serious. From (2.13) and (2.18) it 


follows that 


oP {Q2°ap d | ( l C B 1 S 
O O ‘ O r. A 
Cz | ad (3) i C2 (3); unl or (") ( ») 


The pressure gradient leads to non-zero values of wu and w, leading in 
turn to contributions to (2.12) from the inertia terms which are only 
negligible if 9 


w oes (2.27) 


Hence the viscous terms dominate the inertia terms and so from the 


equation of momentum in the z-direction 


O2q /d ()2 
ws oO\—-| )}. u Ol ae | 
va 


The condition (2.27) then becomes 


‘ , add? pid! 
Q?ad y“B?, i.e. H} he - 


(2.29) 
viotu 
Inserting the experimental values 


H 25,000 gauss 
lor the theory to be applicable. This is much higher than those used in 
the experiments so that even at the highest fields used the inertia terms 
must be greater than the magnetic terms. In the next section, therefore, 
we shall investigate the flow pattern at the other extreme when the imposed 
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magnetic field is zero bearing in mind that in the present section 
c o 

or~ bz 

We shall show that with a similar assumption a complete description 

of the primary flow may be obtained when H = 0, although on theoretical 
grounds it could not be expected to be stable. 


3. No magnetic field 

If a, the radius of the rotating disk in the plane z = 0, were increased 
to infinity, then the flow would be a rigid body rotation with the disk. On 
the other hand, if the radius of the disk is zero the fluid would be at rest 
everywhere. It is assumed, therefore, that for a finite value of a there is 
a region in which the fluid is rotating substantially like a rigid body and 
a region in which it is almost at rest. In between there is a friction layer 
in which the velocity of the fluid changes rapidly. At z = 0 the change 
occurs discontinuously as r increases through the value a and so we assume 
that in this layer, just as in the presence of a strong magnetic field, 


Hence the friction layer is in the neighbourhood of r = a and so in the 
equations of motion we can neglect 1/7 in comparison with ¢/ér. Further, 
from (3.1) and the equation of continuity, uw < w and the equations of 


motion reduce to 12 L ap ’ ‘ 
‘as C c Cc ) 
oe, L—— == @, (3.2) 
a pcr cr C2 
cv ov o-v ow cw C-w l cp — 
e—tw—=—-y—, and w—+w— = y—_—- =. (3.3) 
or Oz cr? cr oz 022 p &z 


The equations of motion (3.2), (3.3) are reminiscent of a boundary layer. 
The main differences are that there is no boundary parallel to the z-direction 
and that the pressure gradient which controls the values of wu and w in 
(3.3) is produced by the layer itself and is zero outside it. This implies that 
the gradients are weaker than in a boundary layer and indeed, as we shall 
see, its thickness is O(v*) as against O(v'). 

In common with the theory of boundary layers the assumption (3.1) 
means that the order of the governing differential equation has been 
lowered so that some of the boundary conditions on v given in section 2 
are redundant and further that the behaviour of the main body of the 
fluid, outside the immediate neighbourhood of r = a, must be known. As 


already explained we assume that outside the layer the fluid is rotating 
almost as if solid ifr < a, whereas if r > a it is substantially at rest. Since 
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> 


the equations (3.2) and (3.3) are wavelike in the z-direction the boundary 


conditions at z = 0 are redundant if w <0 while those at z=d are 
redundant if w > 0. Since the friction layer must be centred on z = 0, 
a we assume that w > 0. The boundary conditions may therefore be 


written as follows: 


whenz=0,w=0 v=Qrifr<a andv=O0ifr>a; 
as Tr > 0.9y> Or. w> 0: 
and as rf > +O, v > VU, WwW YU, (3.4) 


since r does not appear explicitly in (3.2) and (3.3). 
The appropriate solution of (3.2) and (3.3) is found by writing 


u Cys /€z, w = oer, 7 = tn (7 Ri, 
a rd 

m a@a(—| Fs v= QV(n), R= a@Qjrd, 
and P = pa?Q?(z/dR)*P(n), (3.5) 
where F, V, and P are functions of 7 only. They then reduce to 

Cu = : \" F+2yF"}, w an : } F’(n), (3.6) 

z \dR) dR, 
V FV’ = 0, P’ = JV’, 
and F FF" —1F"? = 2(P—nP’) (3.7) 
with boundary conditions 
V(—oo) iy V(co) = 0, F’(+0o0) = 0, P(+0o) =0 
(3.8) 

and P(n)—7 > 0 as » > —o0. 


These equations have been integrated numerically by trial and error. 
First numerical values of F were guessed and the equations for V and P 
ntegrated. The principal difficulty here is that there are two conditions 
nm P so that the assumed form for F must satisfy a complicated integral 

ndition. This was done by adding a suitable constant, found by trial, 
tothe assumed values of F. The third equation of (3.7) was then integrated 

ising the values of P obtained previously and the corrected value of F, 
eading to a new set of values for F, except for an additive constant. The 
procedure was now repeated, the new constant being found from the 
ntegral condition, until the values of F’ obtained in successive trials were 
not significantly different. The numerical work was complicated and many 
repetitions were needed to get agreement even to the second place of 
lecimals. The results are displayed below. 

192.38 L 
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TABLE I 


7 = —- —- .==s = ° I 2 3 4 5 

' an 03 07 15 *31 “48 4 86°59—Cs *56 42 26 14 07 

V IO I'o 98 95 87 “71 “51 *31 17 09 "04 
F(o) — +52, F’” (0) ‘06, P(o) — 63, V’(o) ‘19. 


In Lehnert’s experiments R ~ 2.106 so that the friction layer, of thickness 
about 8 in terms of 7, is about two millimetres thick. 

The approximate solution confirms the assumption made above that 
w > 0. An important boundary condition which is not satisfied is u = 0 
on z = 0; according to our solution w is singular on z = 0. It is possible, 
using Kaplun’s theory (6), to choose a system of coordinates such that this 
singularity is removed; but then both w and w are non-zero on z = 0. On 
the other hand, if we take this induced inflow as a mainstream parallel to 
the plane z = 0a boundary layer can be set up in which @/éz > @/ér and the 
fluid may then be brought to rest at z = 0. The thickness of this boundary 


layer is of order r—a\$(a\12 
af )( ) ‘. 


a d 


and the interaction between the two layers is confined to the immediate 
neighbourhood of r = a, z = 0 where the approximations fail anyway. 

From Table I it is seen that fluid is sucked into the friction layer from 
both sides and is swept up it until it meets the free surface. There it must 
be turned through a right angle and flow parallel to the free surface away 
from r = a in another friction layer. 

This friction layer must also ensure that the boundary conditions on 
u,v, \d won the free surface are also satisfied. From an order of magnitude 
study »f the mass flow and the balance between the inertia and viscous 
terms |): this layer it follows that the thickness of the layer and the fluid 
veloci‘’ in it are of the same magnitudes as before. 

If; > athefriction layer on the surface of the fluid has the characteristics 
of a three-dimensional jet so that it will thicken entraining fluid,until a 
balance is achieved between the flow away from r = a in it and the flow 
towards r = a outside it. 

If r < athe problem is more difficult and the following discussion of the 
secondary flow is speculative. Again the friction layer near the free surface 
forms a jet but this time it is carrying rotating fluid from points of higher 
to points of lower azimuthal velocity. Hence the centrifugal term v?/r in 
the radial component of the momentum equation will retard the fluid in 
the layer, thus thickening it and reducing the mass of fluid in it. The 
variation in v across the layer to effect this is O( R-*), from the radial com- 
ponent of the momentum equation, so that the circulation of the fluid in 
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this layer may well be slow enough to permit the angular velocity to remain 
substantially constant. The friction layer eventually converges on z = d, 
» = 0 where it must turn to flow down the axis of rotation. It is quite 
possible that here the fluid has not quite lost all its azimuthal velocity, in 
which ease a line vortex would occur, but further speculation is unprofitable 


it this stage. 
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The generalization from the present problem to Lehnert’s (1) is again 
straightforward, and we deduce that in this case too the fluid above the 
rotating annulus should be rotating almost as if solid whereas elsewhere it 
should almost be at rest. Lehnert found that the motion was in fact much 
more complicated than this and as pointed out in the introduction, it is 
For H = 8,000 


gauss, however, the magnetic body forces, although still small in com- 


concluded that the present solution is unstable if H 0. 


parison with the inertia terms in the steady equations of motion, ensure 


that the flow is stable in which case the solution is appropriate. 
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TRANSIENT STRESSES IN BEAMS OF 
VARIABLE CHARACTERISTICS 
By A. ROBINSON (University of Toronto) 
[Received 3 May 1956] 


SUMMARY 


R. P. N. Jones (1), considering the case of a beam of constant characteristics, has 
shown that for a disturbance which is produced by a concentrated impulse, the 
bending moment is discontinuous even across the fronts of the shear waves. In the 
present paper this phenomenon is studied for a beam of variable characteristics and 
the law of variation of the discontinuity is determined. The laws of variation for 
discontinuities of the shear force or bending moment across fronts of shear waves and 
bending waves, respectively, which were first obtained by Leonard and Budiansky (2) 
are obtained here incidentally, by a different method. The results are compared and 


interpreted. 


List of symbols 
A _ cross-sectional area 
cg velocity of bending waves (= ,(g/pg) = (£/p)) 
velocity of shear waves (= ,/(7,/p,) = \/(@/p)) 
Dz, density of strain energy in bending (= $8?/7,) 
E  Young’s modulus 
F external force 
G effective shear modulus 
I moment of inertia of cross-section 
k radius of gyration of cross-section (= ,/(1/A)) 
characteristic linear dimension of cross-section 
P_ external impulse 
time 
u _ total deflexion (= wg+u,) 
ug deflexion due to bending 
u,  deflexion due to shear 
x spanwise coordinate of beam 
8 bending moment in beam 
y ratio of squares of bending and shear velocities (= ¢3/c3 
ng stiffness in bending (= F/) | 
HN, stiffness in shear (= GA) 
p Shear force in beam 
o density of material 


(Quart. Journ, Mech. and Applied Math., Vol. X, Pt, 2 (1957) 
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ON TRANSIENT STRESSES 


pg moment of inertia per unit span (= p/) 
p, mass per unit span (= pA) 
4 distributed external force 


w angular deflexion due to bending (= @ug/ézx). 


1. Introduction 

AmonG the papers which have been published in recent years on the 
general motion of a beam, a considerable proportion (e.g. 1-5) is based on 
Timoshenko’s approach (6) in which shear deflexion and rotatory inertia 
are taken into account. It is in fact known that the resulting picture is 


he more acceptable from a physical point of view than the calculations which 
he are based on the classical beam equation (e.g. 7). 


It is to be expected that for a beam which starts from rest under the 


action of concentrated, and possibly impulsive, forces, conditions near 


2 the fronts of the shear and bending waves will be particularly important. 
nd The laws which govern the variation of a discontinuity of either shear force 


or bending moment along the front of shear waves or of bending waves, 
respectively, have been deduced by Leonard and Budiansky (2). They are 
equivalent to equations (3.5) and (4.5) below, and apply equally to beams 
with constant or with variable characteristics. However, in a recent paper 
1), R. P. N. Jones, considering the case of a homogeneous beam (beam 
with constant characteristics), found that for a disturbance which is due 
to a concentrated impulse, the bending moment will be discontinuous even 
icross a shear front. In the present paper we show how this phenomenon 
arises and deduce the law which governs the variation of such a discon- 
tinuity in a beam with variable characteristics (equation (5.5) below). 
Incidentally, we obtain the laws of variation of Leonard and Budiansky 
by a different method. 

Although all these laws are, strictly speaking, derived only for discon- 
tinuities, it is to be expected that they apply also more generally to the 
variation of any other type of concentrated shock which is propagated 


along the beam. 
2. General analysis 


The equations which govern the motion of a beam under the influence 


of a lateral external force, when shear deflexion and rotatory inertia are 


taken into account, are 


2A) 
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o* au ao 
—(ugtu,) = p= = 2%, €), 2.3 
Pome (Upto) = Poe an +O ) (2.3) 
ew ug 
a) <= = —-+ CO, (2.4 
PRO PB Rae ax 


where o(x,t) is the elastic shear force in the beam, f(x,t) the elastic bend- 


ing moment, uw, and wg are the deflexions due to shear and bending, 


o 
u = U,+Ug is the total deflexion, and w = éug/éx is the angular deflexion 
due to bending. Furthermore, ¢(2,¢) is the distributed external force 
and 74, 7g; Po» pg are the stiffness and inertia coefficients due to shear and 
bending respectively (see list of symbols above). The last-mentioned four 
quantities may depend on x but not on f. 


Dividing (2.3) by p, and differentiating with respect to x, we obtain 


Fu, Cw “ 0/1 do é [1 d 
Otea ' at? ~~ aa\p, dx} © ax\p,"]} 


Or, taking into account (2.1), 


1 @o 1 eo l\’do Ow <a]! g 
S ——_ a a a. . a. a —¢ 0, (2.5) 
No A =p, OX" Pp; Ce OF dxt\p, 


( l ) d ( 1 
where it) Sete goes : 
Pol dx a 


Again, taking into account (2.2), we obtain, from (2.4) 


Cw Cw , OW 0 (2.6) 
Pr a4? 7 9 UT >. i - ’ «0 
B ot? Ip Ox" Ip Ox 
, dy 
where Ng = 1B 
dx 


We propose to find out under what conditions certain of the ‘unknown’ 
functions involved in (2.5) and (2.6) may be discontinuous across a curve 


C: t = (x) (2.7) 
in a sense to be made precise presently. For this purpose, we introduce a 
system of coordinates A, » defined by 


A = t—y(2), p= 6 (2.8) 
so that the particular curve C corresponds to A = 0. Then, for any f(z, ¢), 
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) ut J r(x) Fs (yxy 2 
— eax? Mat HO ae 
of _ of , ef ef 


and . 
” a= — OA? ~OAOp * Op? 


Accordingly, equations (2.5), (2.6), when transformed to the independent 
variables A, », become 


| e2 > 72 A2 ] a2 
WwW SS 7 . 


l 
/ (x)) == >t aoa ie —+ 
B®. CA? A? No CACm OACn on, Op” 


| aba | | w (x) (- ) af’ (x  - + of’ (ar) - . $) = Y, (2.9) 
Cu? p Po OA Po 


ind 
nal ey’ (2) 27 £ se L 2oe- Dale PB : a { (ng xb" (x)- gh” (x) ——o 0. 
CA CAC ' Op” = 
(2.10) 
Conversely, we have us(x) u—a, t pe 
e ot 
nd SO us (A S a . is 
Cp Op. 
Hence d Maas | ef at lL @ | of (2.11) 
Ou Ox On Ot Cp b’(a) ox at : 
cand cf l da 
nd similarly . — (2.12) 
; CA us (2) On 


We regard (2.9), (2.10) as a system of linear equations for the second order 
lerivatives €°a/0A*, C?w/cA*. Let A be the determinant of the system, 


l | . ° 
A (oh’(x))* | 
‘ 0 pa—ng(b' (a))? | 
l : . . "= ° 
(us (ar) )* (pg naly (2) ]*) (2.13) 


0, 


nd suppose first that A does not vanish. Solving for é20/@A?, é#w/éA?, we 
then obtain 


{ p Naleys (2x) )*) +2 es. pe 


cA? A\ ’ | No CAC atan Ne Op op? 


| te Cc P 0/1 
- (7) (-} r art te 5(-4)} 


20) C*w * ’ P Cw 
(2pp—— + pp—+ (np (2)-+ mW" (a)) So) 


] 9 2 2 2 A2 
{ - 2 ao , ow l oo A 
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and 
Cw | Cw Cw F Cw 
a —_—__—._ | 2 —+ 9; +. (7 "(a) + now’ (x = o| 
cA? pa— nel’ (x))* PB Neu PB ou? ( ip (2) ip ) CA 


(2.15) 
Now suppose that the properties of the beam, p,, pg, etc., are continuous 
and finite positive, with continuous derivatives across C, and that the same 
applies to o and w and their first and second derivatives, except possibly 
0?a/0A*, C?w/eA2. Then (2.14) and (2.15) show that the last-mentioned 
derivatives also must be continuous across C’. Or, conversely, if the deriva- 
tives €°a/€d? or €?w/eA* possess jumps (ordinary discontinuities) across (, 
while all the other quantities just mentioned are continuous across that 
curve, then we must have 


A [eb’(a ) 17 }( pe nal vs’ (a) }?) 0. (2.16) 
Thus, in that case, either 
] 
(es’ (ar) )? 0, (2.17) 
Ne Po 


yielding a pair of one-parameter families of curves given by 


dx 1 ("*) aa (2.18 
dt ob’ (2) i Bs a 
or pa—ng(b'(x))® = 0, (2.19 


in which case the two families of curves are given by 


= i (22) Lep. (2.20) 
dt us’ (a) A \pg aes 

The quantities c, and cg are the velocities of propagation in shear and in 
bending, respectively. Equations (2.18) and (2.20) define the characteristic 
curves of the problem. 

By integrating (2.14) and (2.15) with respect to A from A etoA=«, 
we can show (see the corresponding procedure for a single equation in 
8, p. 360) that discontinuities of the first derivatives @a/@A, éw/éA as well 
as of the quantities o and w themselves, will occur only across curves 


which satisfy (2.16) provided these functions are (in a certain well-defined 


sense which may be expected to hold in a physical case) limiting cases ol 


continuous solutions. 

In general c, + cg so that the families of curves (2.18) and (2.19) do not 
have any member in common. Hence, at a curve C which is given by 
(2.18), 


PB gl’ (x))? 0, 
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and (2.15) then shows that, provided 0?w/eAéu, 62w/eu?, Gw/CA, and o are 
I I I 


continuous across C, é?w/é\? also will be continuous across C. Moreover, 


integrating (2.15) with respect to A from —e to e we obtain 


Cw 


2 pp L- pp a 1 (ngy” - na’) a) ay. 


| eA | ; | pe—nale )- "PF OO aT cA 


We conclude that, under the conditions indicated above, éw/éA must be 
ontinuous across C’, provided ¢*w/cAépu, 6?w/Cu?, Cw/eA, o remain finite in 
the neighbourhood of C. (A case in which these assumptions break down 
will be considered in § 5 below.) Similarly, integrating (2.15) twice with 
respect to A, we can show that, provided éw/éu remains finite in the 


neighbourhood of a curve C' which belongs to (2.18), w also must be con- 


tinuous across C’. But 
Cw l Cw C7Up wb” (a) Cup $ C*ug 
“ - aa (2) L 
CA (x7) cox if) ca us’ (x) CA cA? 
/ "( r) C{u 
ws ~ w+ ah’ (2) E., 
[ as (a)]* cA? 
. o-u ] CW us” (ax) 
ind sO | w. 


\ us’(ax) CA [ ys’ (ar) |8 


This shows that, subject to the above stated conditions, ¢?u,/éA? also 
will be continuous across a curve C' which belongs to (2.18) and the same 
therefore applies to ug CA and wg and to the tangential derivatives of these, 
“Ile CAC, Cwu Cu". 

We select for further consideration a number of discontinuities whose 
variation is of practical interest. We use square brackets in order to indicate, 
f , 


for a given 1 u, the jump of a function across a given curve C, e.g. 


i 


lim(o, eG). ~<}+ 


€ 0 
3. Discontinuity in shear 


Suppose that the shear force o possesses a jump [o] across a curve C 


which belongs to (2.18 In order to determine the variation of [o] along 


{ 


we go back to the partial differential equation (2.3) and replace o by 
the right-hand side of (2.1). We obtain 


d(x,t) == 0. (3.1) 
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Transforming to the independent variables A, uw (see (2.8)) and pro- 
ceeding as before, we find that 





\ O7u nu ou eu Cu 
_ i 2 4. Se tg a 8 By 
(Po— al" (x))*) —e Pore, Pa Gye | Po aye t Pore, 
By ” = Ou, : 
| pg —P+ (ng b" (x) +004’ (x)) —2— p(a,t) = 0. (3.2) 
op” ran\ 


For a given t = p, write down (3.2) on either side of C, for A = +e, say, 
and subtract. Remembering that 


"ry )2 

Po~ Nol (2 )) 
vanishes on C, while €?wg/@A", ?ug/CACu, and 6?ug/Ou? are continuous across 
C (see the end of § 2, above), we then obtain in the limit 


9 o OU, So ‘g ° 4’ ‘la OUg 0 3 3 
Now, by (2.12) and (2.1), 
Ou, 1 [ou, l [o] 
a ee eee a een 
oA | h(x) | | Ne Ws (ax) 
and so, by (2.11), 
é [éu, C ] l oo] 
AN eae ——5 5}le1— FTL 
Op| ¢ Cu \n, w(x) Nes’ (x) Op 
] 1 . l lc 
~(noV'(x))[o]- = 


: us’(x)(n, b'(x))? da Now (a) Op 
Substituting in (3.3) and bearing in mind that p, =: y,(y%'(x))* we obtain 


7 7.\\ Ole] l d ad ] d cae 
NE tae de IY Mel ay de toh Mle] = 0. 


Multiplying by 7, and rearranging we get 


on, h(a) Ue - he (y,¥'(x))[o] = 0. 
Op Cpu 
But this is equivalent to 
A (oF) 0. (3.4) 
Op \n, b' (a) 
Hence, integrating, and replacing #'(x) by +-c>1, we obtain finally 
“fof? = constant (3.5) 
No 


along the curve C. 
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4, Discontinuity in bending moment 
Next we investigate the variation of a discontinuity in the bending 


moment 8 along a characteristic curve which satisfies (2.20). We have 


CW Cw 
B Na now’ (a — - 4.1] 
(8) = np| | = —npw'@ |Z (4.1) 
Now write down (2.10) on either side of C and subtract. The result is, 
taking account of (2.19), 


OW 


oA 


Cc 
») 
-Pp 


p 
Cp 








(ngb"(x)+ ng’ (x)) Ea = 0. (4.2) 
C 


But, by (4.1), 
C CW l l eB) P 
as aie 4.3 
" | b'(x)( new" ( ax))2 ¢ (np ¥' (7) IB] 1g (x) Cu ( ) 


Substituting (4.1) and (4.3) in (4.2) and sili as in § 2 above, we 
obtain 


2ng h(x) API : (ngv'(x))[B] 0. 
Ou Op 
ind hence = (— aw aU B} ‘ 0, (4.4) 


and finally replacing us’ (ax) by +-Cpa 


Cc 





8|? = constant (4.5) 


long the curve C. 


5. Effect of concentrated external force or impulse 

Suppose that the beam is initially at rest and that at time ¢ = 0 a concen- 
trated external force F begins to act at the point 2 = 0, remaining constant 
thereafter. This case has been considered by R. P. N. Jones (1). Jones 
ses a Fourier transform method and determines the values of the resulting 
Fourier integrals approximately by means of the principle of stationary 
phase. He finds, among other things, that the external force F produces 


nitially a discontinuity in the internal shear force of magnitude [o] = 4F 
neither side of 2 = 0. Now a simple energy consideration shows that the 


nertia forces called into play initially are vanishingly small compared with 
the elastic shear force. It follows that the external force must be balanced 
entirely by the internal shear force and hence that the result just quoted 

1}= 4F) is actually exact (on the basis of the equations of motion 
employed here). Applying (3.5) and taking into account that c, = ./(7,/p,) 
we then obtain for the discontinuity of o along the two curves of (2.18) which 
pass through the origin, 


F 7 7 Go i <4 
a] 5) ("5 ) (5.1) 
la Po 
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In this formula, the upper suffix (°) distinguishes the values of ,, p, at the 
origin. Thus, for constant values of ,, p,, [o] = 4F along the curves 
x+tc,t = 0, which agrees with Jones’s result. 
Let C, be the characteristic curve through the origin, which is given by 
dx ' 
a. 


and let the equation of C, be t = %,(2). Now suppose that the constant 
external force F, which begins to act at time t = 0 as described above, 
ceases at time ¢ = 7. There results another discontinuity, along the charac- 
teristic curve C, which is given by 

t—r us, (2). (5.2) 


More precisely it will be seen that the discontinuity across C, is of the 
same magnitude as the discontinuity across C, but of opposite sign. Thus, 
for given t, as C, is crossed from above, the shear force drops by the amount 
4F (my, Po) 78 p2)* and as C, is crossed in the same direction, the shear force 
rises again by the same amount. For small values of 7 the strip bounded 
by C, and C, is very narrow and so o will vary but little in that region. 
In order to determine the effect of a concentrated external impulse P, we 
let + tend to zero while keeping the product Fr constant at the value P. 
Then C, approaches C, indefinitely, while the magnitude of o in the strip 
between the two curves becomes infinite, after the manner of a Dirac delta 
function. Consider now the time integral of the shear force at any point 
«x > 0. For finite values of 7, the distance of C, from C, in the direction 
of the time axis is 7 throughout. It follows that for small values of 7 the 


t 
integral | o dt changes by the amount 


0 - 
Fr (No Po\t 
2 fa 


on crossing the strip. Thus, in the limit, a concentrated impulse P produces 
t 


a discontinuity —}P(y, p,/7° p®)' in the integral | o dt across C\. 
0 

Moreover, since o becomes infinite on C,, the argument which was used 
in §2 in order to show that the bending moment must be continuous 
across a curve of this type is no longer applicable. To calculate the jump 
[8] we first let + be small positive, and we choose another small positive 
quantity e. Introducing the coordinates A, defined by (2.8) foryds(x) = Yy(" 
we haveA = 0 for the points of C, while, fora givent = , the corresponding 


point of intersection on C, (see (5.2)) is given by 


A = t—7,(x) = +. 
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Wenow integrate (2.10) with respect toA from —e tor+e. Using integration 


by parts in the first term, we obtain 


(pg—naly' (x) ix | | | = (Pp nay’ (x) )*) oy dX + 
+ i. THE 
| (2pg + pg—5+ (ng ¥"(x)- ma (2) =] A | ody = 0. 
, " OAC "© ou* ‘ , CA B 
; —_ (5.3) 


We first let « tend to zero, then +. Then the first and second integrals in 
5.3) tend to zero, but the value of the last integral remains finite (i.e. non- 
vanishing) since o becomes infinite. Thus we have, for small 7, 


4 P 1 
a dn OT (4s a. 
» 


0 0 
\No0 Po 


Also, as « and + tend to zero, the first term in (5.3) becomes equal to 


9 CW CW 
(Pp— NBC, ~) pa( 1 y) : 








CA oA 
where , c2/c2. We therefore have in the limit 
CU Pp 7 s) | ~ 
pa(l—y)|—|+-—= | en 0. (5.4) 
cA 2 \No Po 


But 


Q 


CW 4 l CW Pg vy ] C, 
wale PB “aon y—1)[6), 
| | | | | = 1B Vay’! x ( )[B] 


tore 


und so (5.4) becomes 

“2 (y—1)[B]-4 = (2 Pe} 0) 

c3 2 \n9 8 
a 3} = — 5 (752s) tog. (5.5) 

2 \Ns Po y— I 

Note that according to our definitions A decreases across C' for increasing 
values of x. Hence , 3] (enotes the jump of f as C, is crossed in the direction 
of decreasing x2; for the case of a beam with constant characteristics (5.5) 
again reduces to the formula obtained by Jones. 

In addition to the discontinuities across the curves (2.18) through the 
origin, 8 possesses discontinuities also across the curves (2.20) through the 
origin. The law of variation of [8] along these curves is given by (4.5), but 
our present theory does not furnish the absolute value of [8]. However. 
Jones’s calculations yield 


Pe — 
5} constant = = B (5.6) 
=( 
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for the case of constant beam characteristics, and it can be shown again 
that this relation is exact. But this relation must hold initially for small 
t and x, even if the beam characteristics are variable. Combining this with 
(4.5) we obtain the following formula for the discontinuity of 8 across the 
characteristic curve (2.20), with dx/dt - cg, through the origin, namely, 


0 


[B] nas 5 ("322)* &B - (5.7) 


0 0 0 
No Po : ia l 


In this formula the upper suffix (°) again denotes the values of the 
quantities in question at the origin. Ahead of the characteristic curve, or 
front, under consideration, the beam is undisturbed so that the right-hand 
side of (5.7) is actually equal to the value of the bending moment just 
behind the front. Now 


CW ] 


D; = 38 p2 


Ox 2ng 
is the density of strain energy in bending. Accordingly we may write 
(4.4) as 

c 

(cg Dz) = 9, (9.8) 

Cu 
where D3 is the density of strain energy in bending just behind the front. 
(5.8) is equivalent to 


Ca 
b 


. (cg Dg) ; (cg Dg) 0 
ce. ot 
by (2.11), and this in turn may be written as 


2 Ds ag 8 Dp) = 0. (5.9) 
(5.9) has the form of an equation of conservation for the density of strain 
energy in bending, with cg as velocity of flux. A corresponding analysis 
for the full equations of elasticity (9) yields instead laws of conservation 
for the density of kinetic energy. It is well known that there is a gap in 
our understanding of the connexion between the three-dimensional theory 
of wave propagation in an elastic solid and the beam theory on which the 
present paper is based. Thus, the fact that the velocity of extensional 
disturbances in a beam (which is equal to the velocity of propagation in 
bending, cg) is less than the velocity of longitudinal waves in an elastic 


solid indicates that the presence of the boundaries slows down the bulk of 


a disturbance but no detailed analysis is available of the mechanism by 
which this phenomenon is produced (ef. 10). 





6. Ex: 
Con: 
cross-s 


Then ( 


and sc 
Sinc 


corres 


and S¢ 


Nov 
as [8] 
and h 


the ce 


and s 
cases 
span. ! 
the te 
dent 
as | : 











oan 


‘mall 


witl 


3 the 


rite 











ON TRANSIENT STRESSES 


6. Example 
Consider a beam of homogeneous material and with geometrically similar 
cross-sections whose dimensions are given by a characteristic length! = I(x). 
Then (see list of symbols, p. 148 above) we have 
No GA o [?, Po pA a |? (6.1) 
und so a| . 2 (Pp, Ne)? a 3 (6.2) 
Since |o] is the discontinuity in the total shear force it follows that the 
corresponding shear stresses vary as 1-1, Again 
Ne EAk? o I, 4) »Ak? oc IA (6.3) 
Ip Pp f 
und so, if the variation of [8] obeys (4.5), 


B] & (pg ng)* oc I. (6.4) 
Now according to beam theory, the corresponding normal stresses vary 
l 


is [B\-*, i.e. again as /-!. On the other hand, by (6.1) and (6.2), c, and ¢g, 


ind hence y, are constant along the span of the beam. It follows that for 


1 


the case in which [8] obeys (5.7), 

B| oc (yn, p,)* x l, (6.5) 
und so the corresponding normal stresses vary as /-?. Thus, in all these 
uses, the effect of a variation of the cross-sectional dimensions along the 
span on the stresses is less than it would be under static conditions. Indeed, 
the total bending moment and shear force due to a static load are indepen- 
lent of / and so the corresponding normal stresses and shear stresses vary 
is -3 and /-?, respectively. 
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STRESS SOLUTIONS FOR RECTANGULAR PLATES BY 
CONFORMAL TRANSFORMATION 


By A. M. WINSLOW (University of Washington, Seattle) 
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SUMMARY 

Rectangular plates with distributed forces on the boundaries, x ta, y b, 
are considered. The 2,y-components of the coplanar tractions are expressed as 
monomials or polynomials in positive integral powers of x or y. For these conditions, 
two-dimensional stress analysis given by Muskhelishvili is supplemented by additional 
developments. A function of the tractions, in terms of the complex variable z and its 
conjugate Z, is separated into two parts. In this manner an exact polynomial solution 
is combined with a solution derived by a conformal transformation. The transforma- 
tion, without rotation of coordinate axes, maps z-points of the rectangular area as 
¢-points in a circular area of unit radius. In this transformed solution where finite 
summations are taken as approximate values of infinite series, bounds of the approxi- 
mations can be determined. Typical procedures are shown in examples with square 
and rectangular plates. 


1. Theoretical analysis 
IN using the stress notations due to Muskhelishvili (1, 2), the real and 
imaginary parts of the equations 


Tea TT gy = 2[$3(2) + $(2)], Tre tTyyt ty = 226) (z)+-y(z)] (1) 


7 


furnish a solution for the stress components 7,,, 7 The solution 


rr? yy Try? 
requires the evaluation of two analytic functions ¢,(z) and y,(z) defined by 
the following equations (2) to (5). We introduce Airy’s biharmonic stress 


function U, 


: e2U . ee. 2s (2) 
rr = 7 y? 3 Tyy ox? ; Try = er y “ 
and write 2U = 2¢4,(z)+2¢,(2)+.(z)+¢.(2), (3) 
where W(z) = [ us,(z) dz +-const. (4) 


Equations (2) then show that the constant of integration in (4) can be 
arbitrarily chosen without affecting the derived stresses. 
On the boundary, the tractions satisfy the equation 
s 


,(z)+26;(2)+,(Z) = | (it,—7,) ds +const = T;(z, 2), (5) 


8o 
where traction components 7,,, 7, are positive in the positive x, y-directions 


[Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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ind negative in the opposite directions. In anticlockwise integration around 
the boundary, the initial point s, and the constant of integration can be 
wrbitrarily fixed. The notation 7}(z,Z), introduced by the author, enables 


the traction integral to be converted in terms of x, y by using x = (z+-Z)/2 


and y = (z—2)/27. With polynomial distributions of tractions the function 
[\(z,2) usually has some terms in Z only, which can represent values of #,(2) 


in (5). Also 7\(z,Z) may have other terms representing ¢,(z)+-2¢;(Z). For 


these terms, exact stress solutions for 7,,,, 7 


cx Tyy Try ave Obtained in polynomial 


form either from equations (2), (3), (4), or from (1). 
For the remaining terms of 7}(z,Z), a solution of (5) for ¢,(z) and ¢,(z) is 

lerived by using a transformation 

z= w(C), (6) 
where each point of the rectangular area bounded by x = +a, y = +6 
in the z-plane is mapped conformally on a unit circle in the ¢-plane. The 
wuxiliary complex variable is € = pe*’, and on the circumference of the circle 
we have |C 1. Love (3) develops w(¢) as a power series of the form 

i z Bon aa (7) 

n=0 

In (7) the inconvenient 270-degree rotation produced by the presence of 
the coefficient —7i can be eliminated by changing —i to +1. With this 


modification the transformation, as indicated in Fig. 1, is 


2 oe L— P,(cos 2a)¢2"** | 


where sin a k, bla = K/K’' 


The notation P,, designates the Legendre coefficients, k is the modulus of 
the Jacobian elliptic function snu, and 4K and 2iK’ are the real and 
imaginary periods of the function. Love’s similar transformation to 
8) differs in having K and K’ interchanged and «a replaced by 
x Lir—a, P, (cos 2x’) = (—1)"P,,(cos 2a). (9) 
These changes, combined with the rotation —i, give a result equivalent 
to (8). 
In applying the transformation (6) the following stress notations are 


used. For any function ¢,(z), 


The derivative 


’ a | d(C 
$;(2) d (Cc) — wn ee, (10) 
dz w'(C) 
With o = e*’ representing ¢ on the circumference y of the unit circle, 
092.38 M 
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é =o. Thus a function of ¢ can be expressed as a function of o. The 


following formula, 


l * nd (fn for 7 0 
| 2 do {¢ or n (11) 


271 |} o—C = |0-—s for n < OF’ 
, 

applied to a polynomial 7c, ¢) containing positive and negative integral 
powers of o, gives 


=A0, FDEP _ Tt, t-1)*, (12) 
c—~<¢ 


l { T(c,c—) do 
271 | 

y 
the star indicating that in the polynomial 7'(¢,¢-!)* positive powers only 
are retained. 


With these notations, the boundary equation (5) transformed by (6) is 





$(c) + ale (2) u(e) = T(c,6), (13a) 
@ (Go 
with the conjugate 
Ja) + BOP) 5 60) = M4, 0). (13b) 
w (oc) 


When contour integration is applied to each term of (13a), (13b), 


o(f) + a ale a + const = 7'(Z, f-1)*, (14a) 
271 J w (G)(o—C) 

W(C) + .. ’ mat (c) oF ema = T(f-1,f)*. (14b) 
271 w (a)(o—C) 


. 


In (14a), (14b) the constant terms can be disregarded, since they do not 


affect equations (1). For the contour integral in (14 a), adopting an arrange- 
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ment used by Gray (4), we put 
¢'(¢) y TAY . 
AS) = Bt) = > 6,0". (15) 
w (Cc) n=0 

Differentiating with respect to ¢ each term of (14a) after integration, we 
obtain j 7 
w'(C) B(l) 4 77 (2) BEC 1)}* = —{T(C,¢-1)*}. (16a) 
C — 
In (16a), taking the terms in each power of ( separately, the simultaneous 
equations provide a solution for the coefficients 6, in (15). Substituting 

in (14b) permits evaluation of the contour integral and we have 


(fC) +-{a(lC) B(O)}* = Tl, 0)*. (16b) 
With X(6) = de, 6", (17) 
n=0 


, direct solution for each c, is found from the (”-terms of (16b). These 
values of B(C) and (Cf) substituted in (1), (10), (15) give solutions for 7,,, 


Ty» Tr, a8 functions of ¢ and ¢, 


1] z 
Tre tTyy = 2[ B(C)+ B(C)], | 
sa} ore —, * (18) 
ral se a ee — [o(Z)B'(2)+-H'(0)) | 
w (¢) 
” 2. Theory applied to a square plate 
Typical applications of stress analysis are provided by a square plate 
) bounded by a La,y -a, with tensile forces 7,,. = Sy? on the opposite 
nly th, ' : , ae 
; boundaries 2 -a. In (5), evaluating the boundary integral, we have 
) is iSys 
| (iz,—1,) d. , (19) 
3 
sa 
| T > S >)\3 S 3 5 > 53 
un (z, 2 2—zZ)' - (z3 — 3z2Z+- 322*—23). (20) 
24 24 
sb When (20) is used in (5), three of the terms of 7}(z,Z) give a boundary 
equation g 
$1(2)-+26)(2)-+Hal@) = — 5 (P+ 82829), (21) 
La which is satisfied by boundary values of 
d,(z) Sz3/24, us, (z) S23/24, (22) 
h Substituting the values (22) either in (2), (3), (4), or in (1), we find the 
polynomial solution 
10t Fis Q-] 25 Sax2 + 0-625Sy?”, Tyy 0-375 S2?—0-125Sy?, (23) 


ce -_ - OA Ay 
t Try = U2ZOSITY. 
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These stresses (23) will be combined with stresses derived from the | for7,, a 


remaining term of (20), solution 
ae S2°z followin 

T;(z,z) = —. (24) 
5 strain -€! 


From (6), (12), (24) we have 





mir * S ¥\19 </¥ * a 
T(E, C-*)* = — iLeo(<) Pa(e-")}*. (25) _ 
For the square plate, the transformation (8) is 
2a “ : 
w(C) : — a P, (Q)f4n+1 
‘ 1:854075 <y 4n+1 an(O)S 
n=0 
where for n 0, P,(0) 1, and for n > 0, ‘ (26) 
1.3.5...(2n—1 
Py, (0) = (—1)" nie 
2.4.6...2n 
Using 


= 1-854075 3. Bow 


n=0 n=0 5 
we find the stress solutions summarized in Tables 1, 2, 3. lypic 


# 6. z 5 a : S 2a 2 
Bt) = > 6..c™. Ht) = > cuish*, C= , (27) 








sidering 
TABLE | T 
ABLE In Ts 
; ; The con 
% by Os bye aT boo boy ( 

SS — ee erates steps up 
2°0 —o'l 0°04167 0°02404 001608 —o-'orl172 0°009023 | 1°0270 I I 
1‘ 1°2083 o°1202 0°08042 0°05859 0°04512 0°03612 | —1'0699 For exa 
0°75 0°7163 1°1445 O°1055 OOo 21 0°06501 0°05350 | o-9216 
0625 o*5541 0°0523 I°lI73 0°09390 0°07736 006517 | —o°8186 0: 
05469 O°5117 0°5284 0°6225 I°IOI2 005522 0°07 300 0°7440 _ 2 
cee Z . 4 | po ix-hgui 
0°4922 0°4629 0°4642 0*5000 0°6053 I'0go2 — 007849 0°6871 

o*4512 0°4267 0°4222 —0°4375 0°4824 —0°5934 1°0522 | 06419 each C- 








xccurac 
Table 1 gives the seven equations obtained from (16a) by taking the | an pe , 
“ms in £0 4 £8 f24 senaratelv I 1S , ase equations . 

terms in C°, ¢*, C°...., ¢*4 separately. Values of b,, from these equation: From 

and values of c,,,, from (16b) are given in Table 2. 0-000 


TABLE 2 hree fis 





For ¢ 
n b,,/¢ C4n+g/al 
ia, Ee oe Aaa In (18), 
o | 04847 0° 3942 
I 0°452 0-162 . 
— ‘ f appr 
2 01850 0°092 30 
3 0°1053 006044 Irom th 
| 
06982 0°04242 - 
. 2 a — ries, ¥ 
5 | —0°05065 0°02920 
6 0'03899 progres: 





f the 
Table 3 gives the stresses in the cross-section x = 0. In this tabulation successi 








on 
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for 7, and 7,,,, line (1) shows the polynomial solutions of (23), line (2) the 
slutions of (18). Lines (1) and (2) combined give the total values in the 
following line. Line (e) shows, for comparison, stress values computed from 
strain-energy equations derived by Timoshenko (5). 


«€ 


TABLE 3 





O*4 orf oS I°o 
¢ 430 0°639 0°833 ie) 
2 o'116 255 0°433 
147 ors56 147 0°055 
{o°600| 
7¢ 272 0*402 521 
: lors75 
wi 
7 274 4OS o°521 o°553 
2 SI 087 
5 i] 103 IOI 
) 52 14 fe 


3, Bounds of approximation 

Typical methods of determining bounds of error will be shown by con- 
sidering values of 7,.. for a 0, y 0 and for x 0, y a. 

In Table 1 the numerical coefficients are shown in abbreviated form. 
[he computations were carried out to at least six significant figures in all 
steps up to the final values of Table 3, to reduce cumulative approximations. 
For example, in the third equation of Table 1, 

0-75b) —0°7163466,-+ 1-144761b,—...+0-05355846,, = 0-9216C, (28) 
‘ix-figure precision is readily attained in each coefficient of b,,. In Table 1 
ch C-term, however, involves slowly converging summation and has an 
curacy of +0-0002C. The resulting bounds of error in the stress values 


n be computed by a procedure due to Dwyer (6). 


From (18), (27) for € = 0, 7 2b,. Using Dwyer’s methods we find 
)-:00012C as the range of error in }), giving 7,,, 0-141Sa? correct to 


ree figures for both upper and lower bounds. 


For € = i, wherer,, = 0,7,,,is found by computing the values of r,,.+-7,,, 


6 
18), (23). This requires evaluation of B(Z) > 5,,, having two kinds 
0 


n 

\ipproximation, the error in each b,, and the approximation resulting 
tom the use of only a finite part of an infinite series. Relative to the finite 
ries, we have for B(C) in (27) and Table 2 an alternating series of terms 
togressively decreasing in absolute value, where upper and lower bounds 
ithe summation are found by terminating the summation with two 
successive terms. Moreover, noting in Table 1 the progressive changes in 
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lines and columns of coefficients relative to the values on the major diagonal, 
we deduce the effect of using an additional term in b,, and having 8-term, 
8-equation series. Compared with the 7-term series, we conclude that the 
absolute value of each b,, would be decreased by relatively increasing 
percentages, from b, to b,,, and that the range between upper and lower 
bounds of summation would be reduced and included within the range of 
the 7-term values. This effect of increasing the number of terms was 
confirmed in complete computations with 6-term series (not shown here) 
compared with the present 7-term series. 

Applying the preceding considerations, for b,, we find the range of error 
to be +0-00036C. Noting the range +0-00012C already found for 4,, 
without further computation we take +.0-0003C as an average error for 


Using (18) and Table 2 we have for ¢ = 7: 


6 


Upper bound (7,,+7,,,), = 4| > b,,44 7(0-0003C)| 2 0-1002Sat | 


n=0 (29) 
5 
Lower bound (7,,+-7,,,); = 4| > Den 6(0-00030)| = 0-0753Sa? 
n=0 
Polynomial solution (23), 7,,+-7,,, = 0-500Sa*. (30) 
In the total solution combining (29) and (30), 7,,, = 0, giving 
0-600 Sa2 
T | (31) 


tt ~~ 0-575 Sa2{’ 


as noted in Table 3. 


4. Theory applied to a rectangular plate 

To compare a solution for a rectangular plate with the previous solution, 
we use the same tractions and dimensions except that the boundaries 
y +b are changed by making b = 0-501950a (approximately 0-5a and 
avoiding interpolation in available tables). The polynomial part of the 
solution including equations (19) to (25) is unchanged. For the rectangular 
plate, in the conformal transformation (8) we have 


x 10°. K = 1-582843, K’ = 3-153385, (32) 


»») 


from tables derived by Legendre (7). Values of P,,(cos 20°) up to n = 32 
calculated to six decimal places by Tallqvist (8) were extended to 33, 34 
by an extrapolation formula. Contrasted with conditions in the square 
plate solution, the sequence P,,(cos 20°) is not monotonically decreasing in 
absolute value and is not alternating in sign, providing less convenient 





determ 


B(¢) = 


we obti 
in thes 


Table ¢ 


To 
C-col 
of +( 
gives 
result 
Bi) i 
Using 











STRESS SOLUTIONS FOR RECTANGULAR PLATES 


jeterminations of convergence and bounds of approximation. Using 
\2 
) , (33) 


5, 6. The arrangement 


: 5. ” S 2a 

B(C) > b,, Eo. us(C) y c (an+1 hy = , ree 
= n=0 n=0 8 \3-153385 
we obtain the stress solutions indicated in Tables 4, 
in these tables is similar to the corresponding tables for the square plate. 


Table 6 shows values of T.. i the cross-section x Q). 





b l b,. Cc 








by b, " * a 
I $95 5275 0°02468 0°005531 1°314 
7 . 0°07 404 001659 0°02145 2°521 
+ 9°02700 0°03574 0°07 407 2°553 
37 0°05004 0°1037 0°1296 2°669 
| 
54 0°5007 0°1606 —0°I720 2°140 | 
16 7301 0°73895 OIgQI4 1°295 | 
o°5761 0°7135 08165 05533 — | 
| 
TABLE 5 
be Mt f 
2 “en+1/% 
44 0°4953 
28 02481 
301 Ov11Ig5 
s 0°03803 
1275 O'O00g81 5 
582 2439 
TABLE 6 
\ 6 al £6 
im 482 “686 88 I'o 
3 145 )4 0°400 
Ig "O47 "124 
” JO°375| 
5 341 0°33 
-3| 
7 


To find bounds of error in ),,,, liberal values of approximations in the 
C-column of Table 4 are 0-5533C-+-0-0005C for the last term and a range 
of +0-002C for each of the other C-terms. A solution by previous methods 
gives the average error +0-004C for bo, b, 54,..., by. In applying this 
result to r,,, ¢ i, from (33) and Table 5 we note that the summation 
B(¢) is a convergent alternating summation except for the last term b,, €?°. 


Using this summation up to 6,8 as an upper bound, with the entire 
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summation as a lower bound, computations similar to (29), (30) give, for 


t =i, __ {0:375 Sb? 


Tax = \0-278,Sb2)’ (34) 


as shown in Table 6. Practical considerations indicate that the lower bound 
is not less than the average stress }.Sb* in the cross-section x = 0. 


5. Discussion of results 

In the sequences of Table 5 irregularities in values of b,,, c,, are apparently 
due to the cumulative effect of irregular convergence in the finite series 
(8), (32) for w(f). In Table 2 similar irregular values of b,,, c., were found 
in the square plate solution when only seven terms of the series (26) for 
w(C) were used in preliminary calculations. 

When bounds of 7,.,. in (34) for the rectangular plate are compared with 
(31) for the square plate, it is evident that, for equally small percentages 
of approximation, the solution for the rectangular plate would require the 
finite summation of the infinite series to be taken to more terms. 

In (34) the mean value of 7,., is 0-3265Sb?. This value and the other values 
of 7, in Table 6 differ very little from the average stress 4,Sb*, agreeing with 
results found by Timoshenko (5) in a strain-energy solution. 

Beside reduction of numerical work, an additional advantage gained by 
using a polynomial part of the stress solution is seen by examining the 
values of 7, in both Tables 3 and 6. In line (2), with increasing values 
of ¢/: the percentage errors of approximation are increased by retarding 
the convergence of the power series. This unfavourable trend is reduced 
by line (1) where the exact polynomial solution increases to become even- 
tually a major part of the total solution. 
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THE EQUILIBRIUM EQUATIONS OF THE 
INEXTENSIONAL THEORY FOR THIN FLAT 
PLATES 
By D. G. ASHWELL (University College, Cardiff) 
Received 20 March 1956] 


SUMMARY 
[his paper supplements recent work by Mansfield and Kleeman (1, 2) on the 
ending of inextensible flat plates. The functions of the membrane forces in the 
plates are described, and the existence of shearing forces concentrated along their 
edges is deduced. Equations of equilibrium are obtained, from which the 
magnitudes of the membrane stresses (except near the edges), and of the edge-shears, 
may be estimated. A new form of the differential equation for the deformation of 
such a plate is given in terms of the edge of regression of the developable surface 
nto which it deforms. The accuracy to be expected from the theory, when applied 
practical plates, is discussed. 
1. Introduction 
Ix recent reports, Mansfield (1), and Mansfield and Kleeman (2) have 
investigated the behaviour of thin flat plates that are flexible in bending, 
but inextensible in stretching. The importance of this work lies in the fact 
that, in certain circumstances, thin plates behave approximately as if they 
had these properties—though the degree of accuracy to which this is true 
requires discussion (section 8, below). The inextensional condition requires 
that the plate distorts, under loading, into a developable surface, and the 
letermination of the shape of a given plate under given loads reduces to the 
letermination of the generators of the developable surface into which it is 
listorting. Mansfield and Kleeman have solved this problem by obtaining 
an expression for the strain energy due to bending stored in the plate, in 
terms of the positions of the generators, and maximizing this expression 
give a differential equation which can be solved to give those positions. 
This solution, however, has the disadvantage, common to many solutions 
by energy, that it may not give an adequate picture of the mechanism by 
which the plate supports the load, and may have overlooked features of the 
problem having bearing on design. It seems worth while, therefore, to 
‘supplement the work of Mansfield by seeking a formal solution to the 
problem by a different approach, considering equilibrium and compatibility 


lone. In this paper such a solution is obtained, different in form from 


Mansfield’s, though equivalent to it. 


Before, however, proceeding to the general problem, it will be helpful 


‘0 consider two comparatively simple special cases, each of which demon- 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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strates the existence and function of an action in the plate additional to 
the bending action considered by Mansfield. 


2. Notation 


breadth of rectangular strip. 

curvature ratio, b?/ Rh. 

flexural rigidity of plate, Bh?/12(1—o?). 
Young’s modulus. 

shear per unit length along a generator. 

shear per unit length perpendicular to a generator. 

shearing force. 

plate thickness. 

moment of inertia of cross-section of rectangular strip, bh3/12. 
Yo 
ey : 

| - dy for n = 1, 2, 3. 

Ya ; 

defined by equation (14). 

bending moment per unit length along a generator. 

bending moment per unit length perpendicular to a generator. 
bending moment. 

bending moment given by inextensional theory. 

bending moment given by small deflexion theory. 

lateral load per unit area on plate. 

resultant load on part of plate. 

perpendicular distance from centre of load to generator. 

distance along generator from edge of regression to foot of perpen- 

dicular from centre of load. 

radius of curvature. 

distance along edge of regression. 

tension per unit length along generator. 

tension per unit length perpendicular to a generator. 

angle between generator and fixed direction. 

deviation of rectangular strip from inextensional surface. 

distance along reference axis; a dash denotes differentiation with 

respect to a. 


distance along generator from edge of regression. 


Ya. Y» Values of y at plate boundaries. 


width of region of plate affected by edge deviation. 
, angles between generator and tangents to edges of plate at inter- 


sections with it. 


= \/{3(1—o02)}C. 
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EQUILIBRIUM EQUATIONS FOR THIN FLAT PLATES 
o curvature, 1/R. 

o Poisson’s ratio. 

7 shear per unit length, in the plane of the plate. 


6 P%¢2/DI,. 


3. Pure bending of a rectangular strip 

The simplest problem to which the inextensional theory can be applied 
is that of a rectangular strip bent about the major axis of its cross-section 
by terminal bending moments. In the elementary (small deflexion) theory 
if bending one finds that the cross-section distorts into a circular are, which 
is anticlastic with respect to the principal, longitudinal, curvature, and 
Poisson’s ratio times it in magnitude. The applied bending moment J, 
required for small deflexions is given by 


M, = EI/R, (1) 


where EF is Young’s modulus, / is the relevant moment of inertia of the 
cross-section, and R is the radius of principal curvature. 

The inextensional theory predicts that the cross-section remains flat, 
the strip deforming into a cylindrical surface co-axial with the applied 
moments, and that the applied bending moment required, /,, is given by 

M, EI/R(1—o?), (2) 
where o is Poisson’s ratio. If o = }, this represents an increase of 12-5 per 
cent. over V.. 

Fortunately, an exact large deflexion solution to this problem is available. 
Lamb (3) obtained an expression for the distortion of the cross-section, 
und a more detailed discussion has been given by Ashwell (4). These 
solutions enable some estimate to be made of the condition under which the 
inextensional theory may be regarded as accurate, and this will be discussed 
later; their usefulness at the moment is to draw attention to one of the 
mechanisms by which the shape of a distorted plate approximates to a 
developable surface. 

If the strip were to suffer large deformations of the non-developable 
uticlastie type, predicted by the small deflexion theory, longitudinal 
filaments near its long edges would be appreciably farther from the axis 
of curvature than similar filaments near its centre-line. Thus, in the absence 
of any resultant longitudinal force in the strip, outer filaments would be 
stretched and inner filaments compressed. Since these filaments would have 
i: curvature about the axis of curvature of the strip, the resulting longi- 
tudinal tensions and compressions would have components forcing the outer 
filaments towards the axis of curvature, and the inner ones away from it. 


Thus they would tend to destroy the anticlastic curvature and flatten the 
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strip, as the Brazier effect tends to flatten bent tubes. This effect in plates 
seems to have been first described by Searle (5). In fact, the lateral curva- 
ture of the strip never becomes great, and if the curvature ratio, C (equal 
to b?/ Rh, where 6 is the breadth, and h the plate thickness) is greater than 
about 30, the deviation from the inextensional cylindrical form is effectively 
confined to a region running along each edge of the strip. If 5 is the width 
of this region, it can be shown (Fig. 6, reference (4)) that if o = 4, 


6/6 = 3/VC. (3) 
As C becomes indefinitely large, this deviation, and the resulting longi- 
tudinal stresses, become confined to a region indefinitely close to the edges 
of the strip. 

Since C may be made large by reducing the plate thickness h, and as 
this has the effect of increasing the ratio between the stiffnesses of the plate 
in stretching and bending, it may be expected that a plate for which C is 
large approximates in its behaviour to an inextensible plate. Thus inexten- 
sible plates will resist any tendency to transverse curvature—which may 
be interpreted more generally as curvature about an axis at right angles 
to a generator—by means of stresses in the ‘plane’ of the plate, in regions 
infinitesimally close to the edges. These stresses effectively apply a bending 
moment about the edge of the plate sufficient to destroy the tendency to 
curve transversely. This is one of the two important functions of the 


membrane stresses in a plate which is deforming into a developable surface. 


4. Symmetrical triangular cantilever with tip load—-edge shears 

Next will be considered the behaviour of a cantilevered plate having the 
form of an isosceles triangle A BC, supported along the unequal side BC 
(see Fig. 1). It is loaded by a concentrated lateral 
load P at its tip A. Symmetry suggests that the 
developable surface to which its distorted form 
will approximate has generators all parallel to BC, 
such as EF. Now a surface whose generators are 
all parallel must be a cylinder, and its curvature 
along any one generator is constant. Thus, if the 





cantilever behaves according to the inextensional 
theory, and the deviations from its developable surface necessary to 
prevent anticlastic curvature are confined to regions infinitesimally 
close to the edges, the curvature, and hence the bending moment per unit 


length along EF, must be constant. This bending moment m may easily 
P.AH =m.EF, 


or m !Ptan9@, (4) 


be calculated, since 
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where P is the concentrated load, AH is the perpendicular from A on to 
BF, and @ is the angle ABC or ACB. Therefore m depends only on P 
nd @, and is constant over a given plate with a given load. Since the 
urvature of the plate is proportional to m, the curvature is also constant, 
and the plate bends into a portion of a circular cylinder. 

The condition at the generator EF concerning bending moment having 
been satisfied, the transference of the shearing force P across it must now 
be considered. Since the distortion of a plate into a portion of a circular 
ylinder demands the application of bending moments only, a plate so dis- 
torted can transmit no shearing force; and since the whole of the plate, 
except for regions infinitesimally close to the edges, suffers just such 
listortion, the conclusion is reached that the whole of the shearing force 
is transmitted along the regions infinitely close to the edges. 

This phenomenon, which at first perhaps seems rather surprising, is 
essentially the same as that discussed by Kelvin and Tait (6) when con- 
sidering the behaviour of a rectangular plate acted on perpendicularly by 
1 balancing system of four equal parallel forces applied at its four corners. 
They found that, for small deflexions, the body of the plate suffers dis- 
tortion into equal and opposite curvatures about axes at right angles to 
each other, and that the forces are transferred from corner to corner by 
infinitesimal strips along the edges of the plate. 

In this case, the value of the shearing force F in each infinitesimal strip, 
is clearly, in virtue of symmetry and equilibrium, equal to P/2, but if the 
plate is such that it does not distort into a 
portion of a circular cylinder, such as the //////////////// 


plate in Fig. 2, then it is possible to determine 





the value of F in terms of the local bending 


moment and the angle between the generator 





ind the edge, at the point considered. Thus, 
in Fig. 2, QR is a generator making an angle 
4 with the edge at R, and ST is a second generator infinitesimally close 
toit. 7'U is perpendicular to QR. Let the curvature of the plate as seen 
irom above be convex, and let F be the downward edge-shear which acts 
on the elemental triangle 7UR at T; an equal (sensibly) and opposite 
shear will act upwards on 7'UR at R. Let the moment per unit length 
in the plate about QR at R be m. This is the moment corresponding to 
the curvature at R of the developable surface into which the plate is 
deforming. Then taking moments about UR, for TU R, 


P.7TU = m.U@, 
PF 


m cot @. 
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The derivation of this expression neglects the presence of the membrane 
forces in the plate, but it can be shown that their moment about UR js 
of the second order of small quantities (see section 6, below). Also, since 
the surface is developable, there is no twist in it, and hence no torque on TU, 

The function of the edge-shear F’, then, is to provide the bending moment 
required for the length of generator U R, but it has another effect that must 
be noticed. For by taking moments for 7'U R about 7T'U, it is seen that the 
strip between the two adjacent generators is loaded at its right-hand end 
by a moment of m’, say, per unit length given by 


m’ = Feoté (6) 
or m’ = mcot?é. (7) 


Now the usual anticlastic curvature due to the Poisson lateral strains will 
cause the strip to behave as if it were loaded at its end with a lateral 
moment of om, and thus the presence of the edge-shear F increases its 
tendency to display transverse curvature in the ratio (am-+m’):om, or 
(o+-cot?@):o. This tendency will be resisted in the same manner as that 
in the rectangular strip, that is by membrane forces very close to the edge 
of the plate. 

The symmetrical triangular cantilever has served to demonstrate the 
existence of edge-shears, and the more general problem can now be con- 
sidered. 


5. Edge of regression 

It is proved in books on differential geometry, that a tangent sliding 
along a curve in space generates a developable surface, and that, conversely, 
every developable surface except a cylinder or a cone may be so generated. 
As the tangent slides it coincides in turn with each generator of the surface. 
The curve is known as the edge of regression of the surface; it is the envelope 
of the family of straight lines which are the generators of the surface. 


6. Equilibrium of an element of the plate 

Fig. 3 shows an element of a plate deforming into a developable surface. 
It is bounded on two sides by adjacent generators AB and CD. RR is the 
edge of regression of the surface; AB touches it at O, and CD at O’. The 
other two sides become circular ares, centre O, when the plate is flat. OA 
is y, AB is dy, and the angle AOC is da. m and m’ are bending moments; 
f and f’ are shearing forces perpendicular to the element; ¢, ¢’, and 7 are 
tensions and a shear in the plane of the element, each acting as shown and 
measured per unit length of the edge to which it is applied. p is an applied 


lateral load per unit area. As in the case of the triangular element of Fig. 2, 
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i) 


there are no twisting moments, since the element, being part of a developable 
surface, has no twist. 

The curvature of a developable surface, at a point on a given generator, 
is inversely proportional to its distance from the point of contact of the 





generator and the edge of regression, i.e. in the case, to y. Thus, if 1/R is 
the curvature at A, the curvature at B will be y/R(y+8y). Then, from 
the usual plate theory, 


m’ oD/R 
aDy ? (8) 
{R(y+8y)} 


ind m bm’ 


where D is the flexural rigidity of the plate (equal to Eh3/12(1—o?)). 
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The sides of the element, AC and BD, have lengths equal, respectively, 
to yda and (y+6y) da, to the order of accuracy required, and so, taking 
moments about AC, 

(f'+8f’)(y+y) dady + db py dady?-+-(m’+8m’)(y+dy) da — 
—m'y da —(m+-dm) dyda + by dady?/R = 0, 
From which, using (8) and taking only terms of the lowest order of smallness, 
f'y = m. (9 
The moment about AB of the force ¢+-8t can be shown to be 
(t+-dt)(y da)? dy/2R, 
and thus, taking moments about A B, 

(t+-dt)(y da)” dy/2R+dmdy + 3 p(y da)? dy —(f+4f )y dady = 0, 
giving f 7} (10 
; y\da 

tesolving perpendicularly to the element, we have 
py dady +ty dady/R — dfdy +f’y da —(f’+8f’)(y+sy) da = 0. 
Since, from (9), ydf’+-f’ dy = —ém, this may be reduced to 


t 1 /dm' 1 (df 
pt +2 —= |}, (1] 
R- y\dy y\dx 


Resolving parallel to AC, we have 

dtdy +ry da —(r+87)(y+-dy) da +7 dady = 0, 
dt dr 
——Y 
dx ~ dy 


Finally, resolving parallel to AB, 


giving —27= 0. (12 


(t’ + 8t’)(y+dy) da —t’y da —drdy —tdyda = 0, 
dr dt’ 
— Yy — 
dx dy 


In equations (10), (11), (12), and (13) differentiation with respect to y 


t’tt = 0. (13) 


giving 


must be performed along a generator, and differentiation with respect to a 
must be performed along a line perpendicular to a generator. 

In section 3 we considered one of the two functions of the membrane 
stresses in the inextensional bending of plates. Equation (11) demonstrates 
the existence of a second, namely, the maintenance of equilibrium of 
elements of the plate. The shearing forces f and f’ are called into play by 
the requirement for moment equilibrium of the element, but they cannot in 


general satisfy force equilibrium. Force equilibrium requires the existence 
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f the forces ¢, t’, and 7. Unlike the forces in section 3, they are not con- 
entrated near the edges, but distributed all over the plate. 


7, Equation of the edge of regression 
In order to determine the shape of a developable surface, it is sufficient 


to know its edge of regression, and in this section, the differential equation 





Fic. 4 


for the edge of regression of a plate is derived in terms of its geometry and 
conditions of loading. Fig. 4 shows a typical plate, and two adjacent 
generators which are tangential to the edge of regression RR at O and O’~" 
The angle between them is 5a, measured positively in a clockwise direction. 
The generator OA B cuts the edges of the plate at A and B, where the 
ingles shown, between it and the tangents to the edges, are 6, and 6,; G is 
the centroid of all load on the plate to the left of OAB, and GN is the 
perpendicular from G on to OAB. GN =q and ON =r. s is measured 
long the edge of regression in the direction shown, and OO’ is 8s. 

Since the curvature about a generator at a point on it is inversely 
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proportional to y, and since there is no curvature in a direction along the 
generator, the bending moment at the point may be written 

m = kl/y, (14) 
where k is a constant for a given generator, but varies from generator to 
generator. Then taking moments about A B for the plate to the left of AB, 
we have - 
| om dy = &i,, (15) 


Va 


Pq = 





where P is the resultant total load acting at G, and 
Yo 
Y dy 


I 
y 


(16) 





Va 
The shearing force acting on AB at a point Q is given by (10), where 
the differentiation is along the path QQ’, from one generator to an adjacent 
generator perpendicular to the first. Thus 


(= ld (" 1 /dk k dy 
y\da y dx ,) y*\d ) y? \da) 


But the increase in the value of y in going from Q to Q’ is —6s, so that 


f 1 /dk\  k [ds 17) 
f= a5.) + 5a) | 


Resolving perpendicular to the plate, for the part to the left of AB, 
Yb 
P = | fdy t F - " (18) 


Ya 


where F and F, are the edge-shears at A and B, given, from (5) and (14) by 


a 


a 


PF kcoté, | 


Ya 
elt (19) 
and F, acta. 
Yo 
Thus from (17), (18), and (19), 
P—I (7 kl (= keot@, _ keot 6, (20) 
. d 7 , : d ’ Ya Yn 7 
Uo 1 
where is | oY 
- 2 
- (21) 
Up 
a 
and I, | ey 
« ; y® 
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Finally, taking moments, for the part of the plate to the left of A B, about 
a line through O parallel to GN, 


Ue 


Pr = | yf dy +y,F,+% fF, (22) 
dk 1s 

giving Pr = [,—+kl, + k cot 6,+k cot 4,. (23) 
La aa 


Now k and dk/dx may be eliminated from (15), (20), and (23), and the 
resulting equation solved for ds/da, the radius of curvature of the edge of 
regression at O. Thus, 


ds _— I{4,—rl,|+9|(1,—/y, cot 6, +-(2,—11]yp)cot 4, | 


da ql, I, 1 | 


In the Appendix, this expression is deduced from the differential equation 


(24) 





obtained by Mansfield by an energy method. 


8. The accuracy of the inextensional theory 

Mansfield and Kleeman (2) use the exact solution for a rectangular strip 
in pure bending to indicate the degree of accuracy that can be expected 
from the inextensional theory. This they do by plotting against A? (where 
is 0-41C, if o = 3), (i) dM/dp, where p 1/R, and (ii) the root-mean- 
square value of w when integrated across the strip. In each case they 
draw the curves given by the elementary small-deflexion theory, the exact 
theory and the inextensional theory (Figs. 8 and 9, (2)). The first of these 
two sets of curves is rather misleading, as it suggests that the error of the 
inextensional theory is about | per cent. or less for values of A? greater than 
bout 1-5. But in Fig. 7 of (4), W/p#TJ is plotted against C, for such a 
strip, and the error of the inextensional theory is found to be about 9 per 
cent. when A? is 1-5, about 5 per cent. when A? is 4, and about 2-5 per cent. 
when A? is 13. From these figures, from Fig. 9 of (2), and from the experi- 
mental results given on Figs. 11 and 12 of (2), it is possible to draw the 
conclusion that the inextensional theory will predict the general distortion 
of a cantilever plate of reasonably compact shape, to within an accuracy 
ofa few per cent., provided that the greatest deflexion is at least about 
20) times the plate thickness. 

Other aspects of the problem of the rectangular strip that can be easily 
investigated by means of the exact solution are the area affected by the 
deviation near the edges, the magnitude of that deviation and its effect 
on the stresses. The width 6 of the region affected at each edge is given by 


3), section 3, above. In terms of A this, for a 1. becomes 


5/b 1-92/A. (25) 
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For the typical case of b/h = 100, h/2R = 0-004, giving A? = 33, 6/b becomes 
0-34. Thus only the middle third of the cross-section is effectively free from 
distortion and its associated stresses. If A? = 100, 5/b is 0-19. So in many 
practical cases, the width of the region affected by edge stresses is likely 
to be of the order of one-fifth, at least, of the width of the plate. This 
seems to explain the rather poor agreement between theory and experiment 
near the edges of the plate shown in Fig. 10 of (2). 

The greatest magnitudes of the deviation from the inextensional shape, 
in the edge-regions, are 0-1h at the extreme edge, in a direction away from 
the axis of curvature, and 0-02h just inside the edge, towards it. These 
values cause the greatest tensile stress to be 20 per cent. greater than that 
predicted by the inextensional theory, and the greatest compressive stress 
to be 4 per cent. greater. These figures are for o = }; for o = }, they are 
15 per cent. and 3 per cent. The large error in the tensile stress is very 
localized at the edge of the plate, but it might well be greater in the case 
of plates whose generators meet the edge at angles other than 90°, on 
account of the additional edge-moments applied by the edge-shears 
(section 4, above). 

In the body of the plate, the membrane stresses must be added to those 
given by the inextensional theory. They can be readily estimated in 
particular cases from (11) of section 4, above. 

There seems little doubt that the inextensional theory can be of great use 
in determining the forms of appreciably distorted free-edged plates, and that 
combined with a means of applying a correction for the conditions at the 
edges, should provide adequate data about stresses. This last problem has 
been investigated by Fung and Wittrick (7). 


APPENDIX 
Derivation of equation of edge of regression by energy 
Equation (24) can also be obtained by the energy method used by Mansfield and 
Kleeman (1, 2). They use coordinates « and x to define the positions of generators 
of the plate, where « is the angle between a generator and a fixed straight line, and 
x the distance of the intersection of the fixed line and the generator, from a fixed 
point in the line. See Fig. 5. The distance y, between O and the fixed line can b 
shown to be dz . _ 
Y> sin a. 26 
Y: dx 
The relationship between x and «x to be satisfied by the generators is expressed in 
terms of a function (7, x’, x) (F, in refs. (1), (2)), where 2’ is dxv/da; this, for a uniform 
isotropic plate is 4 Pq?/DI1,. (27) 
The equation obtained by Mansfield and Kleeman is 


-7' - — =~ 


0, (28) 
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itities, regarded as functions of 2, x’, and a. 


dl, 
da’ 


dq 


da® 





(29) 


29) it is necessary first to evaluate, in terms of ¥,, Yp, 
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together with (26) and occasional appeals to the geometry of Fig. 5, it may be 


that 


and 


dl, 
dx 
dl, 
da 
dq 
da 
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shown 


_ sin(@,+a) sin(@,—«) 





Yp sin 6, Yqsin@, 
, cot@, coté 
—I, cos x sin a[cotte , cot Fo), (33) 
Ya Up ; 
‘4 ] ) ; 
= |———Jsina, 
Yb Ya 
- —J, sina, (34) 
sin a, (35) 
0, (36) 
ds 
I, — + cot ™, + cot ,, (37) 
“dx 
ds ecoté coté 
Oy | Micah sw Wott (38) 
dx Ya Yo 
r. (39) 


By these equations, (29) may be readily converted to (24). 
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SUMMARY 


This paper deals with the solution of a mixed boundary-value problem in the 
elementary theory of thin elastic plates. The plates are in the form of half-planes, 
with the boundaries partly clamped and partly simply supported. The solution 
involves the determination of a function, regular in the upper half-plane, the real 
part of the function being given on part of the boundary, and the imaginary part 
on the remainder of the boundary. The method of solution is illustrated by a 


number of examples, some involving isolated loads. 


1. Introduction 
DurING recent years half-plane problems have been considered by Mus- 
khelishvili and Tiffen (1), (2), inter alia, who have given solutions in terms 
of two complex potentials. In determining these potentials, Muskhelishvili 
continued them into that part of the plane not occupied by the material, 
using Cauchy integrals to find them. He confined his solutions to problems 
in the theories of plane strain and generalized plane stress, indicating that 
problems involving clamped or free thin plates were very similar to the first 
and second fundamental problems in these theories. However, problems 
involving simply supported plates have no analogy in the other theories 
and for their solution the method used by Tiffen is preferable (3). 

In this paper a solution is obtained for a boundary which is partly 
clamped and partly simply supported with, or without, the presence of 


isolated loads. 


2. Fundamental equations 
The mid-planes of the plates are chosen as the plane Z = 0 of a rectangu- 
lar cartesian frame O(a, y, Z), and the notation used is that of Stevenson (4). 
The displacement w of the mid-plane in the direction OZ satisfies the 


equation Viw F(x, y), (2.1) 
where F(x, y) is zero if the faces of the plates are stress free. In the case 
F(x,y) = 0 the general solution of (2.1) may be written 

Ww 20)(z)+-20(2Z)+-w(z)+ (2), (2.2) 
where z v+iy and Q(z), w(z) are functions of z which are regular in the 


region R occupied by the plate. 
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The mean stresses £2, yz and the mean stress couples xy, yx, yy (= —: 


tx) 
are given by the equations 
aes sian Oe 
xz+ yz 2P — Vi w, (2.3) 
2 
LY — ya -~P(1+n)Viw, (2.4) 
“9 
~ ~ er > Cw 7 
and xy +yx+2iyy —4P(1—n)—., (2.5) 
ez? 


where 7 is Poisson’s ratio and P is a constant depending on the thickness 
of the plate. 


3. Boundary conditions 


On a horizontally clamped boundary w and @w/én are both zero (5), 
Apart from a possible constant in w, this is equivalent to taking 
Q(z) +-20'(z)+w'(z) = 0 (3.1) 


on the boundary. For a simply supported boundary the conditions to be 


satisfied round the boundary are w = ns = 0. In particular, if a boundary 


consisting of a straight line parallel to y = 0 is simply supported then the 
conditions are w= yx = 0. (3.3 
Using (2.4), (2.5), and (2.1) we find that 


YX 


oP ()’(Z)} 


2(1+-n){Q'(z)4 —(1—n){20"(z)+20"(2)+-0"(z)4 


As in a previous paper (3) it will be found convenient to use the conditions 


(3.4) 
instead of (3.2), where 
l+y Cw 


‘ 2 
> Cx* 


If the boundary is supported under a flexural couple f(x), then the 
conditions to be satisfied on the boundary are w = f(x). 

In the remainder of the paper it is assumed that the plate occupies the 
upper half-plane R, with the line y 
into two parts C, and C,, so that C, consists of the segments a, 
mS 
C, is taken to be horizontally clamped and C, is to be supported under a 
flexural couple f(x). 


- 0 and yz : 


0 as its boundary C; C is divided 
<< b,., 


where k .., nand b, < a;,.,,, while C, consists of the remainder of C. 


The specified boundary conditions are 


ow 


= 0 


Oo <= 


on C,, (3.6) 


cy 
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—f(x 6: oa 
ant Ww 0, v SK on Cy. (3.7) 
2P ' 
\sw — 0on C, 22w/éx? = 0,s0 that the condition (3.7) implies that yx = f(z) 
n (! 


>. 


| 4, A solution of the mixed boundary-value problem 


The plate occupies the upper half-plane R with the conditions (3.6) and 
7) specified along the boundary y = 0. It is assumed that potentials 
(),(z), w)(z) can be found, corresponding to given loading conditions, which 
satisfy w = @w/ey = 0 on y = 0. A general method for determining such 
potentials has been given by Tiffen (2). In general, these potentials will 
give non-zero values of v along y= 0. If these values are given by 
g(x)/2P then, to fulfil the boundary conditions, additional potentials 


\2(2Z), W4(2) are required which satisfy 


Ww 0 on C, 
h(a) (x)+q(x ' 
v “ f(x) +9(@) on (,, (4.1) 
2P 2P 
nd —=0 on.. 
Cy 
This displacement w is expressed in the form 
w 20)(z)4+-202(Z)+-w(z)+a(Z), (4.2) 
cu ‘ jn ’ 
so that 20)'(z)4- 02(2)+-w (Zz) (4.3) 
Cc 
ind v 2" (z) +20" (Z)+-w"(z)+0"(Z) (4.4) 
Consider the potentials ©,(z), w,(z), where 
w,(Z) -202, (2). (4.5) 
‘hese potentials give 
, Cw ’ ” 
w = 0, v frefQ)(a)}, and 4im}{Q,(x)} (4.6) 
cy 


when y = 0. The problem has now been reduced to finding potentials 


onnected by (4.5), and satisfying 


im{Q(2)} = 0 ond, . 
(4.7) 

while re; 02,(x) h(x)/8P onG. 
rhis requires that, on C,, @w/éy 4im{Q,} will be at most a constant. 


The final solution for w must be examined, so that, if necessary, this 


constant can be removed, giving éw/e¢y = 0 on C,. 
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function 
P(z) 


v 


which gives positive 


by means of a Cauchy integral « 
I (2) — 


Hence, if Q(z) is an arbitrary po 
for Q)(z) is given by 


P(z) = {(z—a)(b—z)}'. 


This integral is evaluated by 


whe re 


for a function f(z) 


G,(z), G, 


respectively. 


that This f(z)P(z) 





as 


P(z) is regular in the whole plane, 
values of P(x) as y > 0+ on C,. 
branch of P(z) and reflect it in the real axis, then the resulting function 
behaves as the complex conjugate function of P(z). 


As re{ F(x)} is now specified on the real axis, the function F'( 


anP | 


with poles at z 
infinity, the principal parts of f(z)P 
Re isees 


By using a suitable contour and residue theory, 
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In the determination of the required potentials, we make use of the 


n 


—a;)(b;—2)}*. (4.8) 


and we take the branch 
If we take the other 


cut along C,, 


Thus, defining P(2) 


) in 
this manner, we see that as y > 0, P(x) = P(x) on O, while P (2) — P(x) 
on C, 
The function F(z) is defined as 
F(z) = Q4(z)P(z), (4.9) 
whence F(x)+F(2) {O4(x)+O4(2)} P(x) on C2) 
a 7 . (4.10) 
OF (x)—O}(a)} P(x) on C;) 
Thus, if 4(z) satisfies the conditions (4.7), then 
¢ (h(x)P(x)/8P onC,, 
re{ F(x)} , (4.11) 
; lo on C. 


z) can be found 
f the form used by Tiffen (2). Thus 
ha. (4.12 


C2 


lynomial with real coefficients, a solution 


ea ae } MOPO) gy. 12). (4.13 
87P P(z) . anal P(z) 
The evaluation of integrals obtained from (4.13) 
(a) The integral ' 
| mee (5.1) 


. 
a 


Muskhelishvili (1) for polynomial f(t). If 
the problem is to include isolated forces, this integral has to be evaluated 


Nes, Magers x, in addition to a pole at 


(z) at these poles being 
G,(z), G(2), 
it is found 


L G,(z)+...+G,(z)+G,(2). (5.2) 





(b) 


where } 
The 


| and br 


|= 


t= R 


6. Ap 
tally Cc 
a cons 
flexura 


Hence, 


so that 


where 


42,(2) 1 


where 

The pe 
inspect 
satistie: 
Muskhe 
stamps 


7. The 


The 
tally el 


+. 10) 


t.1) 


ur 


1t101 








A MIXED BOUNDARY-VALUE PROBLEM 


b) The integral 


x 
r f(tyt} . 
J = f( dt, (5.3) 
. 
here f(z) is a function of z regular except at z = a4, X,..., &,, Where tt has poles. 


The principal parts of z'f(z) at these poles are taken to be G(z),..., G,,(z), 
nd by taking the integral along the contour consisting of the circles 


§ and |z R, connected by a cut along the x-axis from x 5 to 
R. we find that 
J /in 2*f(z)+G,(z)+...+G,(z). (5.4) 





6. A plate in the form of a half-plane, with the boundary y = 0 horizon- 








' tally clamped for |x a, and subject to the conditions w = 0, yx = M 
| (a constant) for |x a. The latter condition corresponds to a constant 
lexural couple M. In this case h(a) = M, so that (4.13) becomes 
a 
, iM ’ (a? —#?)! 
(2, (2) —_——__—_——. | = Y dt. (6.1) 
(a2—z?)! 8aP q z—t 
—a 
?)t has a pole only at infinity where 
a” 
a 2 izf] a iz+-O(z-) (6.2) 
Hence, using (5.2) 
* (a°— ; . . » 
it in[{ (a? —z?)*+ iz], (6.3) 
M iM:z iB 
that O%(2) ato (6.4) 
8P  8P(a*—z?)* § (a*—2z?)! 


ere B is a real constant. Unless the stresses are specified at infinity, 


must be o(1) and B 0. Integrating (6.4), 


se Ok 
Q, (a- er, 
SP gP > 
(6.9) 
-™ w,(2 20), (z). 
The potentials (6.5) give im{Q,(x)} 0 on the z-axis if |x| >a. By 


uspection it may be seen that the boundary conditions are completely 





satisfied. A similar result, obtained by a different method, is given in 
Muskhelishvili’s book, as well as by Green (6), for plane problems involving 
stamps, although the boundary conditions are different. 


7, The half-plane under the action of an isolated load 
The plate is to occupy the half-plane y > 0, with the boundary horizon- 
tally clamped for |v! > a and simply supported for |x| <a. An isolated 
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load W is taken to act at the point z Z = a+if8. The boundary condi- 
tions on y = 0 are 

ow 

w= — 0 for |z| > a, 

cy (7.1) 
while Ww yx (= v) 0 for la] <a. 
Potentials corresponding to an isolated load W at z = z are 


w(z) 
where A W /32ahP. 
These potentials may be modified (2) to give w = éw/éy = 0 on y = 0, 
The modified potentials are 


{2,(z) A(z— 2q)log = —, 
ie (7.3) 
W (2) — AZ,(z z»)log — “01 (2)—%,)Az. 
a 
Along y = 0 these potentials give 
- pc 
—— oa (7.4) 





, 2A px ' 12 — t?)8 dt LO (z 
(2) ar a § lt A A = (7.5) 
m(a?—2?)! J {(t—a)?+B2(z—#t) | (a2—z?)! 
a 
The integrand in (7.5) is regular at infinity, but has poles at z = 2p, %), 
where the respective principal parts are 
2 »2)s 2__ 52)4 
G,(z) all ss G,(z) = a =o)" : (7.6) 
21B(z—2») —24B(z — 2) 
The evaluation of the integral in (7.5) by means of (5.2) gives 
Q' (2) 2 Ap iA B(a?—2z?)! 1AB(a?—2z?)! i Q(z) (7.1) 
$24 (2 ones T = 6 i CEL =U T cua 
(2—2)(2—%) | (z—zp)(a®—2?)! (z—2,)(a?— 2)! * (@2— 2?) 
The infinite parts of Q)(z) disappear at z = 29, Z), so that Q4(z) is regular 
in the whole plane, cut on y = 0 from 2 = —a to x = a. In addition, 
(2) (z) is O(1/2?) at infinity, which means that unless the stresses are specified 


at infinity, Q(z) = 0. Integrating (7.7), 


(2, (z) iAB log * zm |! log ((a—z)(a u Zy)}' —{(a T z)(a- “o)}" 
2— iq {((a—z)(a+ 29)}} + {(a+-z)(a—2p)}# 


f > 4 _> \t 
log le a—z)(a+%p)}*—{(a+z)(a —*o)s +C,. (7.8) 








{(a—2)(a+%))}!+{(a+2)(a—2y)}! 











where | 


Thus p 
condit: 


and Q 
The sc 
previo 
tends 


which 
simply 

As 
are w 


functi 


This 1 


and 


In th 
by (7 


8. Ce 

A N 
part « 
soluti 
strair 
at po 
large 
in the 
boun 





ondi- 


x 














A MIXED BOUNDARY-VALUE PROBLEM 189 


where C' is chosen so that Q,(00) = 0. At infinity, Q,(z) = Oif C is given by 


a+2»)'—i(a—z9)* (a+29)'+t(a—Z)! 


(a-+z))!+i(a Z)! (a4 z,)t'—i(a—Z,)! 


(7.9) 


[hus potentials representing an isolated load, subject to the given boundary 
mditions, are (2)(z) 0),(Z) + Q4 (2), w(z) W(z)+-w,(z), where 
W (2) 202, (2), 


ind Q,(z), wo(z), and Q,(z) are given by (7.3) and (7.8, 9) respectively. 
The solution obtained here corresponds to the solution obtained in a 
previous paper (3) for a simply supported half-plane. If the constant a 


tends to infinity, then equation (7.5) becomes 


e 2AR% dt 
: ~ (t—z,)(¢—Z)(z—t) 


. 
x 


(7.10) 


which is the integral giving the additional potentials corresponding to a 
simply supported half-plane. 

A similar solution can be found if the boundary conditions on y = 0 
re w = Cw/cy 0 for x 0, and w » = 0 forz 0. In this case the 
function P(z) is taken to be z! and (7.5) becomes 


0 (2) 2 Arp? _— t* dt \ iB (7.11) 


12 {(t—a)?+ B?}(z—?) 2? 


This integral can be evaluated with the use of (5.4), to give 


Q(z) 2 Ap? _4 ABizy Apizy . (7.12) 
Z—2Zy))(Z—Zpq)  -z2(Z—2%p) -2*#(Z—Zq) 
‘a - ot 1 St o} at) 
and Q), (2) ABi| log ——® 4 log (= + — )|4 eC. (7.13) 
e—é ei — 29 ST 26) | 





In this case Q,(z) tends to C as z tends to infinity, so that Q,(z) is given 


by (7.13) with C (0. 


8. Conclusion 

A solution has been obtained for plates in the form of half-planes with 
part of the boundary clamped and the remainder simply supported. The 
solution is similar in form to solutions obtained by Muskhelishvili for plane 
strain problems by a different method, and contains infinities in the stresses 
at points where the boundary conditions change. The infinities are not as 
large as those obtained in the presence of isolated loads and it does seem that 
in the infinitesimal plane theories infinities may be expected at points where 
boundary conditions change. The solutions obtained by Muskhelishvili 
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and in this paper therefore suggest that at these singular points finite or 
plastic deformations might be expected. 

The author would like to record his thanks to Dr. R. Tiffen, who suggested 
the problem and whose advice has been of great help throughout. 
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THE SCATTERING EFFECT OF A JUNCTION 
BETWEEN TWO CIRCULAR WAVEGUIDES 


By V. M. PAPADOPOULOS 


(Dept. of Engineering, Brown University, Providence R.1.) 
Received 12 April 1956] 


SUMMARY 


The exact solution is found for the electromagnetic problem of the scattering 
effect of a junction between two semi-infinite circular pipes of radius a and b 
respectively (a > b) and having the same axis. The solution involves constants 
which satisfy an infinite set of equations. 

The ease considered is that in which the field is axially symmetric and transverse 
magnetic, and the dimensions are such that dominant mode propagation only is 
possible in each guide at the frequency of excitation. The method leads to a good 
ipproximation to the solution when (a—b)/b is small, and for this restricted range 
the absolute values of the dominant mode scattering coefficients are calculated, and 


numerical results given. 


1. Introduction 

Tue problem to be considered is that of the scattering effect of a junction 
of two perfectly conducting, semi-infinite, cylindrical waveguides of circular 
cross-section, having a common axis; their radii a and b (a > b) are such 
that in the case of dominant EF mode excitation, only the dominant mode 
can be sustained in each guide at the frequency of the excitation. The two 
guides are joined by a perfectly conducting annulus normal to the axis of 
the guides. 

This problem can be discussed in an exact manner either by considering 
the continuity of the field components across the aperture in the plane of 
the junction, or by considering the surface current which must be set up 
on the perfectly conducting annulus in the same plane in order to satisfy 
the conditions of the problem. In each case the form of the incident field 
has to be assumed. 

The first method has been described for the equivalent acoustic problem 
by Miles (1). This method involves the matching of the field components 
ross the aperture, and as a result, an infinite set of equations may be 
found, linear in the scattering coefficients associated with all the scattered 
modes, As Miles points out, this method is satisfactory only when one 
guide has a radius small in comparison with the other; then only a 
small number of these equations need be solved to find a good approxi- 
mation to the solution. In the electromagnetic problem this method 
gives good results when there is no propagation in the narrower guide, 
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but when one or more modes are propagated in each guide the number of 
equations to be solved to reach a reasonable accuracy is inconveniently 
large. 

By using a method which takes account of the discontinuity of the 
magnetic field across the annulus, a result is found which is more amenable 
to approximation, when (a—b)/b is small. The method is general in that 
there may be more than one mode propagated in each guide. The restriction 
to one propagated mode in each guide is one which is used in practical cases. 
It is shown that for small values of (a—b)/b, the current on the annulus, 
which is proportional to the magnetic field discontinuity, and which can 
be expanded in orthogonal series appropriate to the annular region, con- 
tributes to the solution mainly from the first term in the series. The results 
of this paper, in which one junction is considered, may be applied to a set 
of junctions, but only if the distance between each junction is sufficient 
for the non-propagated modes set up at one junction to be of negligible 
amplitude at the next. 

The case to be considered in detail is that with the dominant electric 
mode incident from the narrow guide. The case with the dominant mode 
incident from the broad guide is only slightly different, and this case will 
be described briefly later in this paper. 

The axially symmetric electromagnetic field in this problem satisfies 
Maxwell’s equations, and hence satisfies a two-dimensional wave equation. 
The walls of the guide are perfectly conducting, so that the tangential 
electric field must vanish on these boundaries. The method of solution is 
to assume the incident field and then to use Laplace transforms in finding 
the scattered field. The transform of this field is found on a certain surface 
by the Wiener—Hopf method, and then the total field throughout the guide 
is found, and with it the scattering coefficients corresponding to the incident 
dominant mode. To carry out the analysis it is necessary to assume that 
the free space propagation constant k has a small negative imaginary part 
which later is taken to be zero. The resulting solution can be shown to 
satisfy the conditions of the problem. 

This type of problem has been covered quite fully by a number of authors 
for the junction of two rectangular waveguides. Lewin (2) uses the quasi- 
static method (cf. MacFarlane (3)) to obtain an expression which is then 
improved upon by a variational method. The techniques used in this 
paper, based upon papers by Jones (4, 5) have been applied by Williams 
(6) to the rectangular guide problem. Numerical results in the rectangular 

‘ase are given by Marcuvitz (7) and Williams (6). The only reference 
to the junction of circular guides, that of Miles, contains no numerical 
results. 
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2. The equations of the problem 


The geometric arrangement of the system to be considered is shown in 
Fig. 1. The common axis of the waveguides is taken to be the z-axis, with 


the broader guide, of radius a, lying in the region z > 0, and the narrower 
suide, of radius b, in the region z < 0. Rationalized m.k.s. units are used, 


ind a time dependence exp(iwf) is assumed throughout. Cylindrical polar 
coordinates (7,@,z) are used, but since the solution is axially symmetric, 
there is no dependence on the angular variable 6. In the case of dominant 
B-mode propagation, the only non-zero field components within the guides 








r=a 





are By, E,, E,. These electric field components are given in terms of Bg 
using Maxwell’s equation) by the equations 
1k? _ OB, tk? _ 1 <¢ 
’ and E. = — —(rBp), (1) 
Ww Oz so  * r or 


where k is taken with a small negative imaginary part 


(k k,—tk;, k,,k; > 9), 


o p. i4 a 
1 J2 9 
atts ath — ale = 6. (2) 
Cn” CT Yr or ™ 
Here k, is such that only dominant mode propagation is possible in each 
suide, i.e. 2:405/b < k, 5-520/a (ef. Marecuvitz (7)). This condition is 
only possible ifb<a < 2-295b. 


Since the tangential component of the electric field is zero at a perfectly 
conducting surface, By must satisfy the conditions 


c 
(rBj) =0 onr=a,z> 0, (3) 
cr 
c > 
(7 Bp) QO onr b,z < 90, (4) 
cor 
CBs " 
ind “ 0 onz OQa<r<b. (5) 
C2 


rhe conditions on Bg which ensure the integrability of the surface current 
at the edges of the junction, and which therefore ensure a unique solution 


092 .38 oO 
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(Jones (8), Meixner (9)), are that By = O(1) near the edges r = a, z = 0, 


lé¢ . ; 

and r = 6, z = 0, and —~—(rBy) « FE, = O(c) as the distance o of the 
rer 

point of observation to the edge r = b, z = 0 tends to zero. 


The scattered field is in the form of an outward travelling wave with 
respect to the origin. With k; > 0, this field satisfies the radiation con- 
dition. Apart from the direction and magnitude of the incident field, the 
problem is now specified completely and uniquely. 


3. Notation 

The following notation is used. The zeros of the functions J,(«a) and 
J,(xb) determine a series of real positive constants «,,,, and «,,,, in increasing 
order of magnitude for all integer n 1, so that 

I (K na @) J (K np b) 0. 

Similarly the positive constants «,, in order of magnitude are defined for 

all n > 1 by the zeros of the function J,(Ka)¥)(«b)—Y,(xa)J,(«b), so that 
J, (x, a)Yo(«,, d) Yj(K, AJq(«,, d). 

We note here that as a > b, (a—b)x,, > nz for all n > 1 (10). The asso- 


P,»», and p, are defined by the equations, for n > 1, 


k?)}, 


ciated constants p 


= 
Pna = (k2,—k?), Pny = (x?2,—k*)!, and p, = (x? 
where that branch of the square roots is taken which gives these constants 
positive real parts. The restriction on k, in section | is such that im(p,,) 
and im(p,,) are O(1) as k; > 0, while re(p,,,) and re(p,,) O(k;). The other 
constants defined here have real parts which are O(1) as k; > 0, and 

imaginary parts which are O(/;). 
It is convenient also to define the constant p, = ik, and then the corre- 


sponding constant «x, is zero. For large values of n, 


7 7 1) 7 l 
Ky cle O ; Knog = —(n 1) . | ; Knp (n— 1) ! ol ; 
: a—b n ‘ a n b n 
nar | 7 L' 
so that Pn Lo ) p. (n —})+ O ), 
a—b n a n 
and ee 7(n 1) 410 
. Pnb h 4 n ’ 
The functions ¢,,(7) are defined, for all n > 1, by the equations 
, (7) ae Ji(ky r)Yo(K a)—Vy(ky Joke a), 
] 
and for n = 0 by ¢,(r) = -.- 
, 
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\s defined for n > 0 the functions ¢,,(r) form a complete orthogonal set in 


‘, r b with weight function r?. 


4, Solution for the dominant electric mode incident in the narrow 





éuide 
The first case to be considered is that with the dominant mode incident 
he . rm: ° . > 
nthe narrow guide. This assumed mode has its magnetic field component 
in the form J,(«,),, 7)exp(—p,, z), Which is the form of the dominant mode 
t up in an infinite guide of circular cross-section and of radius 6 (cf. 7). 
\ total magnetic field By which satisfies the conditions enumerated above 
taken in the form 
"8 Bo(r, z) Hr, z)+4,(«y,r)exp(—py,z), in r < b for all z, (6) 
nd B,(r, z) “Wr,z) ina>r> 6b forz > 0. (7) 
| Both 4(r,z) and ‘3(r, z) satisfy equation (2). Jn addition, from (1) and 
1) it follows that é(r4),.-; 0 for z < 0, and from (1) and (3) it follows 
| that @(r“¥3),_,,/er = 0. The functions ‘Y, and ‘Y, are exponentially damped 
| away from the origin, with a damping factor of order k;. It is therefore 
SSO ssible to define Laplace transforms of these functions, so that 
yb4(8,? | ‘W(r, z)exp(—sz) dz, 
D nd yb,(s | ‘V(r, z)exp(—sz) dz, 
ther ae, 
, ere $18 a complex Variable. 
ind . . 
From the equations satisfied by Y we see that 
C* l < a l : 
rre = of ke a}vals, r) 0 forr < 6b. (3) 
cr r or rm 


‘imilarly from the equations satisfied by ‘Y,, and from the boundary con- 


| ) dition (5), we see that 


Sa ae Wo(8,7r) = 8(Bo).-9 = sF(r), say. (9) 
Or? ron 2 
Here x = (s?+-k®)!, where that branch of the square root is taken for which 
0 when re(s) > 0, and F(r) is a funetion proportional to the current 
the annulus z O,a ? b. 
The edge condition (B, O(1) as the point of observation approaches 
0,r = a, or2 0, 7 = b) ensures that 7} F(r) is an integrable function 


ir,and therefore that r! F(r) may be expanded in a series of the orthogonal 
lunctions r'¢,,(r). From the continuity of By with respect to r follows the 


ontinuity of the transform of B, corresponding to the region z > 0 across 
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the surface r = b. It also follows from the integrable behaviour of @(r By) /ar 
near r = b, z = 0, and its continuity elsewhere, that 0(ry,)/ér and @(ru,)/éy 
are continuous at r = b. 

Since the solution of (8) which is bounded at r = 0 is given by 

w,(s,7) = const J,(xr) 
,-p/Or, we can see that 
, 
(reb,)p Jy (Kr) 
KbJ,(«b) 


in r < b, and writing (ry,), for ry,(s,r)| 


ub, (8,1) 





(10) 


Similarly the solution of (9) for a-> r > 6 for which (r,)), = 0 is given 
in terms of (r.);, by the equation 
[ Jo(Kb)¥y(Ka) —Yo(xb)Jy(xa) |bo(s, 1) 


r 


: Asm Jo(Kb)¥, (Ka) Y, (Kb) Jo(xa) )| [ pF (p){J,(«p)Y,(«r) -Y,(«p)J,(xr)\ dp js 
b 


a 


+-48a| J,(Kr)¥o(Kb) —Y¥, (xr) Jo(xd) | ( pF(p){J,(kp)¥ (xa) —Y,(xp)Jy(xa)' dp+ 


b 
i (Ea) (er (wa) -¥i(«r)Jq(«a) |. (11) 


For r = 6 equation (11) becomes 
[ Jo(xb)¥o(Ka) —Yo(Kb)Jy(xa) jybo(s, b) 


‘i ead (cb) Yala) -Y,(«b)Jq(«a) |+- 


pF (p){J,(Kp)¥o(xa)—Yi(kp)J,(xa)} dp. (12) 





Kb, 

b 
Now let , = ( ‘Y(5, z)exp(—sz) dz, 

0 

. 
and p= | V(b, z)exp(—sz) dz, 
so that us, Ly = us, (s, b). (13) 
Similarly let (rus...) = | (7'¥;), exp(—sz) dz, 

4 

and (rf_)’ = | (7r'¥,), exp(—sz) dz, 
so that (ryb,.)’ + (rp_)’ = (rdy)h 


Then the functions 4%, and (r,)’, %.(s,b), and (nba are functions of 8 
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regular in the half-plane re(s) > —re(p,,,), and the functions _ and (ry_)’ 
sre regular in the half-plane re(s) < re(p,,). 

From the continuity of the total field and its normal derivatives across 


the surface r b, z 0, we have that 
Jk 
i} 1(Kyp 5) bo(8, b), (14) 
ST P1p 
nd (wy = (rebo)). (15) 


From the boundary condition (4) it follows that 
(rys_)’ 0. (16) 


From these conditions and equations (12) and (10), we see that 


Lic. by | PE (adel (<p) Yolo) —Ya(xp)Io(wa)} dp 
J3\K 4,9) b 


(rds G(s) : — — —<———— — AS (17) 
st+py Kb J, (Kb) ¥,(xa) — Yq(«b)J,(«a) 


; 2 J,(Ka 
where (7(8s) _ ) 


7K*b? J, («b)| J,(Kb)¥,(«a) ¥,(xb)Jo(xa) : 


(18) 


The funetion G(s) can be written in the form CG, (s)/G_(s), where C is a 


mstant and 


l+s 
s+ik 4 \(1+-s/p,,)(1+s Pn) 
nd where y(s) : ‘alna—blnb—(a—b)In(a—b)}. 
) is a function regular in the half-plane re(s) > —re(p,,), and G_(s) 


1 function regular in the half-plane re(s) < re(p,,). The functions G(s) 
nd G_(s) are derived in Appendix I, where it is shown that, as |s| > 00 in 
he half-plane re(s) re(p,,,), G(s) ~ const(s~!), and that 

y (s)G ( 8) l. 
Returning to equation (17), multiplying this throughout by G_(s) and 
rearranging the terms, we have that 
Altea te 
lh G (s 1\ Ky, 9) GQ (s G_{( Pw») 
S+ py 
rd, )'CG,(s) Sy ( yp 6)G Pw)| 
8 Pio } 
where 
| pF'(p)| J,(Kp)¥o(xa) Y,(xp)Jg(xa)] dp | 
P(s) G_(s)\— 





Kb Jy(Kb)¥y(Ka) Yo(xb)Jy(xa) i J 


the left-hand side of (19) is composed of two terms, one regular in the 
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half-plane re(s) < re(p,,), and the other regular in re(s) > —re(p,,). To 
split the right-hand side into two such terms we proceed as follows. From 
Cauchy’s theorem, ‘ 
: l Ptr 
P(s) - PS) ag. (20) 
271 ¢—s 


where s is any point within C, and where C is a finite rectangle with two 
of its sides parallel to the imaginary axis contained within the strip 
re(Piq) < re(<) < re(pyp). 


Now as |C| > © within this strip of regularity, P(¢) = O(¢-!). The inte- 


grand in (20) is therefore O(¢-'), so that if the rectangle C is taken to be of 


infinite length, then the contribution to the integral from the ends of the 
rectangle is vanishingly small, and hence P(s) may be expressed as the sum 
of two line integrals each of which is absolutely convergent. Thus, 


Cot+i« Cy+i@ 
I , Ake oe l P tt) sp 
P(s) ; G) dt - ~— dt, 
2m J ¢C—s8 271 J  C—s 
C2—ix C\-—ia« 
where -re(p,,) < C, < re(C) < C, < re(p,,). 


The first integral is regular in re(s) < re(p,),), and the second in 


re(s) > —re(pj,). 
Equation (19) may now be written in the form 
Cetin 
e L [i + AD. 
ly (s)G_(s)— Aer a (e)—G_(—p,)]|— — | (5) ge 
| 8+Dyp, 2mtu_ J Ss | 
Cz—in 
Cyi+i« 
J,(«,,5)G i 3 
6 | oG.(s)(rb,)' 22k» )) ( Pi) 1 | P(o) acl, (21) 
| stpy 2m J t—38 
C,—i« 


where the left-hand side is regular in re(s) < re(p,,), and the right-hand 
side is regular in re(s) > —re(p,,). There is a common strip of regularity, 
so that each side of (21) is the analytic continuation of the other. Each side 
is therefore equal to an integral function. 

Now the left-hand side of (21) is O(s-!) as |s| > 00 in the left-hand half- 
plane, since %_, the transform of a bounded function, is O(1/s), and the 
integral { {P(£)/({—s)} df = O(1/s). The right-hand side is similarly 
O(1/s) as |s| > co in the right-hand half-plane, since the Laplace trans- 
form of the edge condition shows that 


(rb,_)’ ~ s-*! = const 


as |s|—>0o in re(s) > —re(p,,). This is because as z>+0+ on r=), 





~~ 








(rBa)’ 


theref 


The i 


Il 


um 


poles 


(rbs,.)’ 


where 
and 
It nov 
that 
u,(s8 
(8, 


Similk 


It car 
simpl 
can a 








A JUNCTION BETWEEN TWO CIRCULAR WAVEGUIDES 199 





To rBo)’ O(z*4). The integral function to which each side is equal must 
‘on therefore be zero. In particular, we find that 
(29 led. 1 tf PO gr, 1 Ale ®)C(—Pw) (99) 
(Vy CG (Ss) Dri 3 ¢—s = CG (s) s- Pw 
‘WO | The integrand in (22) has simple poles in the half-plane re({) < C, at 
py and at Z p, for alln > 0. Evaluating the residues at these 
poles we see that 
’ a 
nte ‘te ] |G_( 4p )dy(K b) = Xn r 7 
vl; ) a Pip)¥1\"109) | S = pl'(p)¢,,(p) dp + 
e of UG ,(s)| 8+) “4 S++ Pp, 
die 1 i 
the a 
Xo y | or 
un + - | F(p) dp}, 23) 
87 Po. } 
where Yo {b In(a/b)G( py) 


TK y Iy(Ky a) 
n 2hG \P, VI y(k, b)| tI y(K, a) R 


0 


ind X 





. 21° 
(x, 5)! —1} 
It now follows from (10), together with (23), (16), and the equation 


yt (ref_)’ (rib,)). 


that 
J,(xr) (G anda (K ype = Ya 
ta oh L ) | ( Py (Ky, 9) L S Xn pF (p)¢,,(p) dp\. (24) 
KbCG.(s)| $+ Diy, <4 8+Pn . } 
b 
» Similarly, from (11), (23), and (15) it follows that 
)[ Jo(xb) ¥o(xa) —Y(xb)J,(«a) | 

(2] bsr| Jy(Kb)¥q(xa) — Yo(xb)Jy(xa) | pF (p)iJ,(xp)¥, (xr) —Yy(kp)A(«r)} dp+ 


b 


} s7r| J, (K? \Y,(«b) Y, (Kr) Jq(«b) | . 


rl eg 

= pF(p) J, (kp)¥o(«Ka) Y(xp)Jy(xa)} dp =. 

- l | J, (xr)¥,(«a) Yi («r)Jo(«a) | 

the UG. (s) Kb ; 

arly } G_(—DiJ, (1p, ~ i 

| | Peet > | pF (p)bn(o) dp\. (25) 
‘ans | > Pib - s- Pn | 


It can easily be seen that u,(s,7) is a function whose only singularities are 
simple poles at s 1D», 8 Pna for all n > 1 and at s = —p,,. It 


} 


can also be shown that y,(s,r) has for singularities only simple poles at 
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8 = —p,,foralln > 1. We are now left with the problem of determining 
the function F(r). Let us expand r*F(r) in orthogonal series so that 


Then for all integer n > 0 


An Ny — [ pF(p)¢,(p) dp, 





b 
a 
where N, = | plbn(p) |? dp. 
b 
Thus N, = In(a/b), 
2 - \\2 
and Ny a Si. (aie) ! for n > 1. 
m2 \ Jy\ Knp) 
It follows from (23) that 
(rb aa : (G ( bs oe Kib pb) an Xn Nn). (26) 
CO)\st+py | Ly stp, | 
Now since #,(s,6) is a function regular in the half-plane re(s) > —re(p,,) 
so that the function «*xh(s, b)[ Jg(Kb)¥,(«a) —Y,(«b)J,(«a)] vanishes at all the 
points s = p, for n > 0, it follows from (12) that for s + Dn 
(ra), = — | pF (p)d,(p) dp 
2/0 ,,(b) . n 
Pn ay Ny (27) 
$,(b) ‘ 


The continuity condition (15) enables us to set up from (26) and (27) an 
infinite set of equations linear in the coefficients a,._ These can be written 
for all nm > 0 in the form 
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This set of equations determines the unknown constants a,. The inverse 
transform ol (8, r) is defined by 


Cot+ia« 
; 
7 » ae hener ay > 
H(r, 2) = ys,(s, rexp sz ds, 
re 1x 
| where re(Py_) < Cy < re(py,). 


This integral can be evaluated by the method of residues if we close the 
ine contour of integration with a large semicircular are contained in the 
region re(s) C, or re(s) > C, according as z > 0 or z < 0. On the semi- 
ircles the integrand is then vanishingly small and its contribution to the 
ntegral vanishes. Thus for r < b, from (24) and the expansion of F(r) 
when 2 0 

<,(7;2) > RY Ti (Kmp )EXP(+ Pinp Z)s (29) 


where for all m 
RG | (G i Prv)y(K i 8) LS tame N, \ (30) 


" > ' a 2 
( BD npy(Kmp 0)G (Pino) Pio T Pmb Pn T Pm! 


n=0 
nd when 0 
. : = . > 71) ‘ Co = ‘ 
Hr, i expl Pry 2)A (hy 7) > rc Ji(Kna T exp(— Pma2)s (31) 
mn 1 
where for all m 
m4 bJy(k maO)Km ! |G ( -Pyp)J (yp 5) S Ay Ly N,,) (32) 
2 2/7 a I : bide 
a Pma Ji Kma a ) | G Pm | Pip Pma Boge Pn Pma 


From (25) and the expansion of F(r) the inverse transform of (sr) may 
be worked out in the same way. It is found that fora > r > b, whenz > 0, 


J(Kina rjexp( Pma)- 


n 


Y(r,z) = > TS 
l 


The expressions on the right-hand side of (29) and (31) represent the 


scattered field for a dominant electric mode incident in the narrow guide. 
5. Results for the dominant electric mode incident in the broad 
guide 
In the case with the dominant mode incident from the broad guide, we 


issume an incident wave in the region r < a to be of the form 
EXP( Pq 2)Iy(Kia r3 


Therefore, by taking scattered fields, represented as before by the 
lunctions ‘Y\(r,z) and ‘Y4(r,z), the total field is 


Bg = Vy(r,z)+exp(p,,2)4,(Kiar) inr < 6b for all z, 


V(7,z)+exp(py.2)44(k,,7) Ina >r > bforz > 0. 
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The solution in this case is found as before. The function F(r) is expressed 
in orthogonal series form 


F(r) In Prlt), 


n 0 ~ 
and the infinite set of equations corresponding to (28) which is satisfied 
by the constants b,,, is given by 


x 


p bin(Sinn T Cont : E,; (33) 
m=0 
where E - Pn(O)K rn bJo( 1. ab) (Pra) 


Pr(Pna Pra) (py)N 
The similarity between the equations (28) and (33) may be noted. The 


expressions which result for the scattered field are of the form 


s T' J (Kmnpy T)CXP(+ Pp 2) When z < 0, 


Kmb 
m=1 
and > R® K(kmaT)eCXP(—Pmaz) When z > 0. 
m=1 
Here 
T\?) ae ; (CK 14 DI y( 445) )G (Pra) he 3 b,, Xn N,, \ (34) 
‘ b 2 Dmp* J, (x Kmb b)¢ 'G, (Pp) | Pimb Pia soot. Pro Pmo) 
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R®) —— le ma b)K ma = >! b, a Xn 4 Nn ; C¢ Rig a OI y( Ky qb)G, (Pi a)| 
Pina @ i(k ma @)}°G_( : Pina) n Pie Piat Pma 
(35) 


6. Approximate values for the dominant electric mode scattering 
coefficients 
A brief discussion of the numerical solution of the problem is now possible. 
The computation can be divided into three stages. The first is to evaluate 
the infinite products which arise, both in the infinite set of equations, and 
in the expression for the scattered field. The second is to solve the infinite 
set of equations, and the third is to sum the infinite series which represent 


the scattered field amplitudes. Unless only one or two specific results are | 


required, an electronic computer is required to do the work. No estimate 


of the error involved will be attempted here. 

When ¢« (= (a—b)/b) > 0, the problem is greatly simplified. We know 
from Appendix II that G,(p,,), G.(py,), and G,(p)) are O(1) and that 
G.(p,) is O(e) for n > 1 ase>0. Also p, = O(1/e) and x, = O(1/e) for , 
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in (28) can be written explicitly in the form 


‘ . 12 
Jy Kma NIo(K mb )j 
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m? 1L® ( 2 - <a aq 
( BD (K mnp) Ji\ Kaplt- (Jo(Kna)}? lem Pr\Pr T Pm)GA(Pm) 
rm, n 0 and as « > 0 this expression is O(e?). When m 0, 
( Iy(* ma) 


3 Oe). 
When 0 
Uy 
hy O(1). 
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* OU 22 
we—>U., 


nly a) is non-zero, and all 
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O(1) for all n. 
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If in the infinite set of equations (28) we assume a solution such that 


the coefficients a,, = 0 for m > 1, then from 


0) we have that 





the first equation of the set (n 
A(1+-Co9) Dy. 
Jf The (r+1)th equation (n can be written in the form 
a, D, Ay ( or = Am mr? 
? m=1 
(" s0 that applying the first approximation of the solution to the right-hand 
(35 side we have that - D,—ay Cy, 0(1). 
> 
. Now we can find a second approximation to a, from the first equation of 
“ing iar 
6 the set, and this is 0 
Y 
j Dy > D.C r0 
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ible 0 0 
1a ( 00 > ( ro Or 
' r 
wn D 
ana - ! Ole). 
nl 00 
St If we now consider the scattering coefficients 7) and R, we can see that 
are ' by neglecting all but the first term in the infinite series for F(r), that is, 
na by using the first approximation to the solution of the equations (26) we 
ire neglecting all terms of order greater than the first in e. A similar result 
nov holds for the equations (33) where we have a solution 
th E, 
b, = Ole), 
{ 0 ] O 
} ~“00 
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Thus if we neglect terms of order « in the infinite set of equations (26) and 
(33), we can find the dominant mode scattering coefficients correct to order ¢, 

On applying the first approximation to the scattering coefficients and 
putting p, = ik throughout, we have that 


RY — A vh- Pu) x 
: 2b*p?, C[G,(Py)P\ik+ py) 





{22a B OAL. EYP Gh Po Gb—DyP 
= 1—2C ae k*[ G(tk) |? .’ 


Ri?) = * [dol Kia b) aa Pia)! x 


2a*p? | Ji(K1,4) |? ik—PDr4 























, {1—2C In(a/b)b?k* >| (tk) Pl (tk—p,_)/(tk+-p,,)}*) 
~ | 1—2C In(a/b)b?k>{ G, (ik) |? .’ 
ma) — bIy(K 1a) Kia G4 (Pra) — Pra) (tk—pyp) e 
, BD qIy(K iy 9)(Piy—Pra)F(Pry) (th—pyq)(th+- py) 
, (L=2Cn(a/b 7G. (id) (it Pradik+Du)/(ik-+Pra)ik—py)) 
*| 1—2C In(a/b)b?k?| G, (tk) |’ 
and TY Pia T. (36) 
Pp 0” 
To evaluate the behaviour of the absolute values of these scattering 
coefficients, as « > 0, we need the results of Appendix II that 
l 
G(s) - A. 
(8) ik (€) 


for finite values of s in the right-hand half-plane, and that 
G k? k , 
F+(Pra)| _ 4 4 (2 , | Oe), 
G.(Pyp) Pi” | Paw), 


together with the formulae 
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If we now put x,,/k = Q, where under the conditions of the problem 
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[@ is proportional to the wavelength A of the excitation, and is given by 
Q = 0-346A/b.| These results are obviously valid when « = 0. It should 
also be noted that (36) is a statement of the Reciprocity Theorem so that 
these final expressions for 7} and 7’ satisfy this theorem. Graphs are 
given in Figs. 2 and 3 to show the behaviour of the first-order terms in the 


scattering coefficients. 


7. Method of solution for other modes 

In this paper we have solved completely only the case with one axially 
symmetric 7M mode capable of propagation in each guide. As long as we 
consider only axially symmetric modes, either 7'E or 7M, there is no 
difficulty in considering the cases with other modes propagating. It is 
necessary in solving such cases to take each incident mode alone and to 
find the scattering coefficients associated with it. The results can then be 
superimposed, 

No work has yet been done in the general case when the field has angular 
dependence. However the variable @ can be eliminated by a finite Fourier 
transform, so that we have a modified two dimension problem which can 
be solved by the same methods as in this paper. 


Acknowledgement 


The author wishes to acknowledge the helpful suggestions received from 


Mr. E. Wild and from Professor M. J. Lighthill during the preparation of 


this paper. 


APPENDIX I 
The function G(s) defined in (18) can be expanded in infinite product form in the 
following manner. The function is an even integral function of s which can b 
expressed in terms of its zeros by the equation 


J,(Ka) J)(ka) | | (1 5, expt) | | (1 = 
1 


x 
va’ 
1 


sa 4 
Jexp{——) (cf. (11). 


The choice of the exponents in this expansion is not unique; any exponent is suitable 
for which the separate products converge absolutely, but the final solution of the 
problem is independent of this choice. Similarly we can write 
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ie | function G can now be written as the ratio of infinite products, in the form 
‘ e)/G s). wh rr 
sae ( 2 J, (ka) 
uf mb? J,(kb)| Jy(kb)¥,(ka) — Yo(kb)Jy(ka)]’ 
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It follows in particular that L,(s) is bounded as |s| + ©0 in the half-plane of regularity 
of the function G,(s). Using Stirling’s formula we find that as |s| > 
G(s) 


Pte) ~ const|s texp|x(s)— [alna blInb—(a—b)In(a »)}\|. 


If then we put x(s) s[alna—bInb—(a—b)In(a—b)]/7, we ensure that 


G(s) ~ consts i 


as |s| > oo in the half-plane re(s) > —re(p,,), and from this definition of y(s) we can 
see by inspection that G,(s)G_(s) 1. It also follows that G_(s) ~ const st 
s| > o in the half-plane re(s) < re(p,)). 


as 


APPENDIX II 


The products G,(piq), G(pip), and G,(p,) for all m > 0 are implicit not only in 
the formulae for the dominant scattered field, but also in the infinite sets of equations 
(26) and (33) which define the coefficients a, and 6,. In this appendix we are con- 
cerned with the behaviour of these products as e€ (= (a—b)/b) > 0. 

The function G,(s) is bounded for any finite value of s in the half-plane 


re(s) > re(Piq) 
when «€ > 0. It can also be shown that G,(s) 1/(s+ik)+O(e) as e > O for the 
same range of values of s. Thus 
; I . l 
Gi(Pra) — + O(e), G (Py) ——_—_— +1. (}(€) 
Pratth Py t+tk 
I 
and . G (ik nm he Ole). 
(7k) ik (e€) 
Let us consider the form of G,(s) at the points s = p, for all » 1, where for 


small values of ¢«, p,, no/eb+-O(e). We know that Z,(s) is bounded at thes 
points so that |(p,+7k)G,(p,)| behaves like 
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This follows from the equation (A2). Thus for small values of €, using Stirling’s 


formula we have that 
\(P1,+tk)G,(p,)| ~ const n! exp[n(Inn— 1)] 


~™ const nt, 


for large values of n. Therefore G,(p,) O(e) for any n > lase >0. 
Finally, we need to evaluate the product |@,(p,,_)/@,(p,p)| as € > 0 and to do this 
we start from equation (A 1). If at this stage we put k; = 0 then we may write 
Pia i|Pyq| and Pib i|Prp|- 


We now have that 
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A NOTE ON THE DIFFRACTION OF A DIPOLE 
FIELD BY A HALF-PLANE 


By W. E. WILLIAMS (The English Electric Co., Ltd., 


Guided Weapons Division, Luton Airport) 
[Received 17 May 1956] 


SUMMARY 
A simple compact solution is given for the diffraction of the field of a dipole by 


a perfectly conducting half-plane. 


1. Introduction 
THE solution for the diffraction of the field of a dipole by a semi-infinite, 
perfectly conducting, half-plane was first given by Senior (1), who solved 
the problem by the integration of a plane wave solution over a suitable 
contour. Two more papers have appeared since on this topic, one by 
Vandakurov (2) and the other by A. E. Heins (3). 

The method of solution that immediately comes to mind is the differen- 





tiation of a suitable scalar solution; closer examination, however, shows 
this approach is not feasible since it introduces too high a singularity in 
the field components at the edge of the half-plane. Vandakurov avoids 
this difficulty by writing each field component as the sum of two terms, 
one of which is regular at the edge whilst the other is singular there. On 
applying conventional transform theory he obtains a solution for the 


yo ~Ts—“>>—— 


regular parts of the magnetic field components normal to the edge and 
also for the total’magnetic field parallel to the edge. The condition of 
vanishing divergence is then applied to calculate the singular parts of the 
field components. Three different and perpendicular directions of the 
dipole axis are considered, and for the axis parallel to the edge of the half- 


— 


yo 


plane a vector potential solution is used instead of the approach described 
above. 


: , ait 
Heins attacks the problem in a different manner; he sets up integral 


~ 


representations for the field components and by comparing these integral 
representations with those for the scalar problem he succeeds in writing 
some of the field components in terms of the solutions to the scalar problem } 


(cf. 4) and an unknown function. He then obtains a Poisson’s equation for 


this function. Heins also considers three perpendicular dipole directions 

separately. ' 
The object of the present note is to present a simple and compact method 

of solving this problem which involves very little manipulation and which 


4 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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enables a general dipole direction to be treated directly. The method used 
ssimilar in some respects to that of Vandakurov but is more compact and 
requires no consideration of three different dipole directions. 


2. Formulation of the problem 
The half-plane is taken to occupy the region y = 0, x > 0, all z; for 
convenience we also define polar coordinates (r,@) by x = rcos6,y = rsin@. 
The dipole is taken to be at the point (29, yp, Z)) or, in polar coordinates, 
», 09,29), and its axis is taken to have direction cosines (1,m,n). If we 


yssume a time variation e’” then we have the equations 

curl E iw, H, div H = 0, 

curl H iwe, E | i mj 1 nk |d(x Xy)d(y- Yo)0(Z—Zp)- 
Here Ee’, He’ are the electric and magnetic fields (in m.k.s. units) 
respectively, and €», x» are the dielectric susceptibility and magnetic per- 
meability in vacuo; 8(a) is the Dirac delta function and i, j, k are unit 
vectors in the x, y, z directions respectively. 
From the above equations we obtain 
(V2+k?)H eurl| /i mj nk |d(x Xy)O(Y Yo)0(z- -Z) 

curl, /i mj nk |d(x ry )O(Y Yo)O(Z —2Zo); (1) 
where k2 = we). and the subscript zero refers to differentiation with 
respect to (2. Yo, 2)). If a new vector A is defined by H = curl, A, then A 
Ils seen tO satisfy 

(V2+k7)A | i mj nk 


The boundary conditions for H on the half-plane are 


O(x X_)d(y Yo)O(z Z)- (2) 





H. O8, 
i, = — P= @ 
, cy cy 
In order to solve this problem we need to define two independent 
solutions d. us of 
(V2+k?)u = 8(a—2y)(y—Yg)8(2— 2) 
| 
. CW 
such that on the half-plane ¢ 0. 
cy 
If we now take A, = 1d, A, = ms, A, = nd then equation (2) is satisfied 


— . Cc 
together with the boundary condition on H,. If we add —n — (¢—y) to H,, 
; CYo 


( 
j 


(¢—u) to H., then equation (1) is still satisfied and the boundary con- 


ditions on the other magnetic field components are now satisfied. The only 
condition of the problem not satisfied is divH = 0, so we have to add 


suitable functions to the field components to satisfy this. These added 
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functions must be solutions of the homogeneous wave equation, vanish or 
have vanishing derivative on the half-plane, and satisfy the radiation con- 
dition at infinity. Thus for these solutions to be non-zero they must be 
unbounded at the edge. The formula obtained for H, represents this func- 
tion as a combination of the derivatives of % with respect to x, and y, and 
hence H, behaves in a similar manner to ¢ at the edge; and since 4%, being 
the solution of the scalar problem, is bounded at the edge, so is H,. Thus, 
since one of the boundary conditions on H, is that it is bounded at the edge, 
we need only consider adding functions h,, h, to H,, H, respectively. 
We thus have to determine two functions h,, h, 


,» Such that they satisfy 


the equation (V2-+k2)u = 0 


and éh,/éy = h, = 0 on the half-plane. They have also to satisfy the 


subsidiary condition 


chy , thy _ lo Ob _ @ eb) é (e+ ©) — a (4 }. ; 


+m 
éx by \ex Cy, Cy Cay} Gz \Ox Oi Gzy\Cy  CYy 


From the known solutions ¢, % we have 
~H®(kt){ /r r j 
Pe 4 23-2 (kt) 2.4 leas d0c0s 36, f. 
OL OXy t \WWr ro) a = j 


op | cd a H}?(kt){ ’o Zs loin 14 sin 10 q 
oe 8% co wer’ rol : west 


Op eb HP (kt) /r ir , 
4. =- 41 (RO) [Fo \cos 30 sin 30, = h, 
Cy CY, t \|Nr COC ' 

where ¢? = (r+r,)?+(z—z,)? and H{?(x) is the Hankel function of the 
second kind of order unity. It is easily verified that f, g, A all satisfy the 
source free wave equation. 

If we seth, =AD+h®,h, = h®+h and define hY = —nh, ho = ng 
ah  ch®). ;, ; 
—*_4 4 jis equal to the sum of the second and third terms on the 
6x oy 


then 


right-hand side of equation (3). This gives 
2}, (2) 5}, (2) ‘ . . (2)/]. 
oh’ ch - ae 7 ; H‘?)(kt 
—#_ 4.9. — 4- 24 len cos 46, —l sin 46,} — cos 46. 
oy oy Czg\N 7 ro) ss P t 
This equation was obtained by Vandakurov; its solution is easily 
obtained and is given by 


’ i : 
hY aa (m cos $6)—Il sin $6,) — F 0. (kt) cos $4, 
\ (7 ? 0 : “0 
” i : e ; 
Ae) — (m cos $6,—I sin $6) —- Hj) (kt)sin 30. 
= "a 
V\TT9 “0 
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Inspection of the formulae for h®, h®, h®, h® show that they are of 
order r-$ at the edge; the other parts of H,, H, are derivatives of ¢, 
with respect to 2, Yo, 29 and so are bounded there. It then follows that 


the formulae obtained here for H,, H,, H, satisfy the required conditions 


it the edge and therefore represent the unique solution of the problem. 
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DIFFUSION FROM AN INFINITE LINE SOURCE 
LYING PERPENDICULAR TO THE MEAN WIND 
VELOCITY OF A TURBULENT FLOW 
By T. 8S. WALTERS (University College of Swansea) 
[Received 6 October 1954.—Revise received 7 June 1956] 


SUMMARY 
A solution is obtained of the problem of diffusion into a turbulent atmosphere from 
an infinite line source of constant strength which lies at right angles to the mean wind 
velocity, the source being at height h above the ground. The mean wind \ elocity U(z 
and the coefficient of diffusion K,. in the vertical direction at height z are taken to be 
Uz) = Ue / end K, = Az. 


Roberts (see, for example, 1, p. 73), in an unpublished paper, obtained the solution 
for a source at ground level; here the general case is considered of an elevated source. 
Roberts’s solution being derived from it as a special case. 

Expressions are obtained for the distance 2), down-wind (from the elevated source) 
to the point at which the ground-level concentration attains its maximum value Xu 
and for the value of yy. 


1. Introduction 

Let an infinite, cross-wind, line source of constant strength Q be at height / 
above the ground, in a turbulent atmosphere in which the mean wind 
velocity is a function only of height z above the ground, the direction of 
the mean wind velocity everywhere being at all times parallel to a fixed 
direction. The source is supposed to emit, at a constant rate, matter in 
the form of fine particles whose presence does not affect the flow properties 
of the atmosphere, @ being the mass emitted per unit time per unit length 
of the source. We choose an axis Ox, at ground level, parallel to the 
direction of the mean wind velocity, and an axis Oz vertically upward. 
Denoting by U(z) the mean wind velocity at height z, and by K, the 
coefficient of diffusion in the vertical direction, the differential equation 
for the mean concentration y(a,z) at a point (x,z) in the two-dimensional 
case in the steady state is 


ne ae Of 5k 
U(z) 4 (x! " (1.1) 
Cx O2z\ co 
The expression we assume for the mean wind velocity is 
U7, 2m 
"(z ait 1.2) 
Ue) = +, 
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where U, is the velocity at some reference height z,, and m is an index 
shose value lies between 0 and 1. We make the further assumption that 
fusion from the elevated source is described by the coefficient of diffusion 


K, = Az. (1.3) 
The value of h must be such that the diffusion zone lies within the region 
f validity of the power law representations (1.2) and (1.3). On the usual 
sumption that the horizontal Reynolds eddy shearing stress 7 is inde- 
pendent of height it may be shown (1, p. 72) that the value of the index « 
1.3) is | 


whose value may be calculated from a knowledge of the meteorological 


m. - Azl-™, where A is a factor 


The resulting expression K, 
factors in any particular case, is now generally accepted as an empirical 


rule in the interpretation of field trials. The calculations are here carried 
ut without restriction on the value of « in the first instance, and the 
solution when a |—m follows as a special case. The problem is solved 
n section 2, the method being, briefly, as follows. 

Let the source of the entity which is being diffused be the infinite strip, 
perpendicular to the wind velocity, whose section by the plane Ozz is 


0, h h 
The presence of the strip is not intended to have any effect 


a, the concentration at all points of this area being 

(constant 
nm the flow, and is to be regarded as merely a mathematical device to 
obtain the result as a0. The solution for the mean concentration at 
points down-wind of this source is obtained in the form of an infinite integral ; 
ind the solution for an elevated line source of strength Q (expressible in, 
for instance, g. sec.~! em.~} units) at height A, whose length is perpendicular 
to the mean wind velocity, is obtained from this by making a > 0 and 
) > © in such a way that the product ay, remains finite and equal to a 
certain constant times Y. As a special case, the solution of the problem 
of diffusion from a line source at ground level is obtained, and is the same 
is that quoted without proof by Roberts in an unpublished paper (1, p. 73). 

In section 3 the ground level concentration down-wind of the elevated 
source is considered, and the maximum ground-level concentration, and 


the point at which this maximum value is reached, are found. 


2. The infinite cross-wind line source 


The transport equation (1.1) becomes, on inserting the values of U(z) 
ind K, given by (1.2) and (1.3), 
1 é, mM 


oo. (-=2x) 2 
Azy’ ox Oz\ 
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The boundary conditions are 


(i) x/xo>1 whenx>0,h<z<h-+a, 

(ii) y>O whenza>0,0<2z<handz>h-a, 
(iii) y>O whena>oorz>o, 
(iv) z"(éy/éz) +0 whenz>0,2> 0. 


(Condition (iv) is the no-flux condition, which expresses the fact that the 
ground z = 0 is impervious to the entity being diffused.) The function y 
also satisfies the down-stream flux condition 


| Uy dz = Q, the integral being independent of x. 
0 


We now make the substitutions 


. Azm #2 
§ = —i—-z, {= 29, (2.2) 
U; 
where 6 = 14+3(m—a). (2.3) 
Equation (2.1) then becomes 
CX ex _ 28+ 1 ex (2.4) 
ae C be ra a ; 
where 8 1(a—1)/6. (2.5) 


The equation is now expressed in terms of the Hankel transform 
RE) = [ x (po) dl, (2.6) 
0 
where J.(x) denotes a Bessel function of the first kind of order s. Proceeding 
in the usual way, by multiplying equation (2.4) throughout by the kernel 
(s+! J. (pt) and integrating with respect to ¢ from £ = 0 to £ = o, we obtain 


the differential equation m 
dx 


dé 


the solution of which is x = O(p)e-$", (2.8) 


px, (2.7) 


where O(p) is a function of p which has to be determined from given 
conditions. 
In (€,¢) coordinates the boundary conditions at € = 0 are 
/. {Xo in h,<€<h,+a,, 
x™\o inod< C<h, and (€>h,+4,, 
as € > 0, where h, = h? andh,+a, = (h+-a)*. The condition to be satisfied 
by x(€) as € > 0 is readily found from the above, and is 


x> lth +4), 4(phy + pay) —hi**J,,,(phy)]. 
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Hence the expression for X is 


x th -ay)8*1,, (ph, + pa,)—h§*1J,, (ph) je$?". (2.9) 


The value of y is obtained from ¥ by using the inversion formula 


x 


" pxdSp) dp, 





x= ¢* 
0 
which gives 


, xo 6-* | [(hy+4y)8 (hy pa, p)—hit (hy p)J,(Cp)e$”* dp. 

0 (2.10) 
This is the expression for the concentration at points down-wind of an 
infinite strip whose plane is vertical and whose edges are at heights h and 
h+a above the ground, the direction of the mean wind velocity being 
perpendicular to the plane of the rectangle. We now calculate the value 
of the flow Q per unit time, per unit length measured parallel to the edges 
of the rectangle, from this source. It is 


} 1 h+a 

, ; uy. £ ro U, 

Q ¥9 U(z) dz Xo 1 zm dz Xo “1 [(h+a)™+t—hm+1], 
; ie a ; (m+ 1)zz 
m+1)z™¢ 
1.€. [(h a) Amttly. = ~ aw. a (2.11) 
1 
Equation (2.10) may thus be written 
(m L)z"’ @ = 
X U. th sand jmtl) [(hy+-a,)**1d,, (hy p+a, p)— 


l Je 
0 
he+1J.. (hy p)|J(Cp)e$”* dp. (2.12) 
The solution for a line source follows from (2.12) by taking the limit as a, 
tends to zero It is 


x 


" cd Ep? > 15 

X = Tr paps | Pela P)A(Sp)e*”" dp. (2.13) 
L> 4 
( 


The ground-level concentration in the case of an elevated line source is 
now readily obtained from (2.13) by allowing ¢ to tend to zero. It is 


g 
x 
bz" @ 
28I"(s- L) Uy 5 


0 


x(x, 0) pst d(h, pe Ep" dp, 


which, on using a result (2, p. 393) for the infinite integral, becomes, finally, 


In 


terms of the original coordinate x and the height A in (x, z) coordinates, 


x(a, 0) = Ba-Be-ve, (2.14) 
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where 
: . -oyn({2m—at+s' 
B= (m+1)QU,/(m—a+ 2)? ABgmM4-Dim—a+ 2] (eos) 
M—A-+- So 
m+1 1 tates 
B ; 9? v= 9\2 4+m’° (2.15) 
m—a+2 (m—a+2)?Az? 


The result contained in equations (2.14) and (2.15), and some interesting 
deductions from it, are considered in section 3. 
For a source at ground level the solution may be found by taking the 
> e ‘ z 
limit as h, > 0 in (2.13). The result, in terms of the original coordinates, is 
; QU e-Uita- x + 2) exp| U, zm—Aa+2 {(m— x-+ 2)2Azm™ at] 
: (m eal 2)m +ax)/(m—Xx-+2) 4 Pzma—1)(m—a+2) 1B) a8 


-— 
4 











(2.16) 


which is the result obtained by Roberts. The result in the case « = 1—m, 
the one usually taken by meteorologists, may be obtained from (2.16). 

An expression for the concentration down-wind of an elevated source 
may be found from (2.13), on using a result (2, p. 395) for the infinite 
integral in (2.13). Expressed in terms of the original coordinates the result 
finally obtained is 


Qha—wl2 ais n | ; ~ miioaiins | Dry [ z\(m—x+2)/2 

(2,2) = — 1 t-witexp| —7%/14 (2 eft | 
xX, 2) A(m—a-+ 2)" exp| x\ h Hl B iF h 
(2.17) 


where J,(x) is the modified Bessel function of the first kind of order v. 
We now verify that this solution satisfies the down-stream flux condition 


( U(z)x dz = constant = Q. 
0 


Denote by F the value of the flux. Substituting from equation (2.17), and 


making the substitution f — x(m—a42)/2 (2.18) 
we have : 
F = 2CQhA-wi2e-Ch™= : | 1B, CPT, _[2Ch™ x+22¢] dt, (2.19) 
0 
where C = yh-*1*-*z-1, (2.20) 


Using a well-known result (2, p. 394) for the infinite integral in (2.19), the 
result F = Q is finally obtained. 

Loberts’s result (2.16) for the concentration down-wind of a line source 
at ground level also follows from (2.17) by taking the limit as h > 0. 


3. The ground-level concentration down-wind of the elevated line 
source 


The ground-level concentration down-wind of an elevated infinite cross- 
wind line source is given by equation (2.14). Denoting by x, the distance 
measured along the ground to points where the ground-level concentration 
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ittains its maximum value y,,, it follows immediately from (2.14) and 


215) that 


iu 


Ty Km-—a+2 
U, hn 





Inserting this value of 2,, in (2.14) gives 


(m 


x-+ 2) 


° 3.1 
(m+1)(m—a- 2)Az ( ) 
ee (3.2) 





ePI}(2m—at-3)/(m—a+2)}U, hm 


The ground-level concentration increases from the value zero at points 


vertically below the source to its maximum value y,, at points distant x, 
yway, and then falls off to the value zero at infinite distance away. The 


value of 2,, varies as h” 


1. For the average value of m in adiabatic conditions 


pecoOMes Ly, x h2 


I 


2. If the value of « is taken to be 1—m this 


n the atmosphere, which is approximately }, this result becomes x, oc h}*. 


The value of y,, varies inversely as h”+! (it is not without interest to note 


that the index a does not appear in this result); when m = } this becomes 


,c ht, Also, x,, varies inversely as U}. 


We now write the expressions for 2,, and x,, in the case when a = 1—m. 


The empirical expression used for A, in computing atmospheric diffusion is 


K. a(n)v" U} n . n(1—n)/(2—n) (2—2n)(2—n) | (3.3) 
n aerodynamically smooth flow, where n 2m/(m+1). In (3.3) v is the 


kinematic viscosity and a(n) is a parameter given by 


an) = — 


9 arl2n\l—n 
(2—n)(zk?n) (3.4) 





(1—n)3 2n93 3n? 


ere k is von Karman’s constant. In rough flow, Sutton replaces v by N, 


1 


the ‘macroviscosity’. The 


value of a(n) for k = 0-4 and for values of n 


between 0-10 and 0-49 is given by Sutton in an appendix to his Atmospheric 


Turbulence (1, p. 104). Identifying our factor A with the expression 


ain yy" U} 


1e result for a y becomes 


ey 


n 4+—n(1—n)/(2—n) 
~] 


9 27 7n B(2+n)/\(2—n) 
(2 n) l th 


\ n2/(2—n) 


2(2+-n)a(n)v"2z4 


(3.5) 


where, to conform with the above, we have worked with n rather than 


with m. The expression for y,, when a = 1—m = (2—2n)/(2—n) is 
Xu C(n) QU; ‘@ (2—n) f 2/(2 n) (3.6) 


where O(n) 


1.0.G. Sx >r'TON, 


4.G. N. Watson 


(2—n)(2+-n)?/2@+Me2 


{1 tmosphe ric 


A 


Tre atise 


9(4+n)/(2+n) 
2) 22+) Ps (3.7) 
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TWO METHODS FOR THE NUMERICAL SOLUTION OF 
MOVING-BOUNDARY PROBLEMS IN DIFFUSION ANp | 
HEAT FLOW 


By J. CRANK (Courtaulds’ Research Laboratory) 


[Received 26 July 1956] j 


SUMMARY 


Two methods are described for the numerical treatment of heat-flow problems in 
which a transformation boundary moves through the medium. There is a corre- 
sponding problem in diffusion through a medium containing fixed sites on which 
some diffusing substance is instantaneously and permanently immobilized. In the 
first method the problem is transformed by a change of variable from one involving 
the simple heat-conduction equation with an awkward moving boundary to an 
eigenvalue problem in which the equation is slightly more complicated but the 
boundary is fixed. In the second method, Lagrangian interpolation formulae are 
used to develop finite-difference approximations to space derivatives of temperature 
based on values at points unequally spaced in the direction of heat flow. This enables 
the course of the transformation boundary to be followed. 


1. Introduction 
WE consider a problem of diffusion in a medium containing a number of 
fixed sites on which the diffusing molecules can be permanently immobilized 
and withdrawn from the diffusion process. If the withdrawal takes place 
instantaneously the problem is characterized by the presence of a sharp 
boundary surface separating a region in which all the sites are occupied 
from one in which none are occupied. The concentration of molecules free 
to diffuse is always zero in the region of unoccupied sites. The boundary 
surface moves through the medium as the diffusion process supplies the 
molecules needed to occupy fully successive layers of sites. 

There is the corresponding problem of heat flow in a medium in which 
a latent heat of transformation is absorbed or liberated at a definite tem- 
perature. Thus, if heat is supplied to the upper surface of a block of ice 
in a container, a boundary surface, on which the temperature is 0° C., 
moves down having water above and ice below. The latent heat corre- 
sponds, of course, in the diffusion problem, to the capacity of the attractive 
sites for the diffusing molecules. 

Numerous attempts have been made to deal with this mathematical 
problem in particular systems. Carslaw and Jaeger (1) give a formal 
solution for the case of a semi-infinite medium of which the surface con- 


[Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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entration or temperature is constant. This type of solution has since been 
veneralized by Danckwerts (2). Stefan (3) and A. V. Hill (4) pointed out 
that in certain cases an approximate solution is provided by the steady- 
state concentration or temperature distribution which would be obtained 
it each instant if the boundary were fixed. Pekeris and Slichter (5) have 
evaluated the first-order correction to the steady-state solution for the 
problem of ice formation around a long cylinder whose surface temperature 
is variable but always less than 0° C. Eyres, Hartree, and others (6) applied 
finite-differences to follow the course of the moving boundary. They con- 
sidered heat flow in a plane sheet which was thought of as comprising a 
number of layers. The temperature of the layer in which the transformation 
is taking place is assumed to remain constant at the transformation tem- 
perature till the whole layer has passed through the transformation. The 
method has since been used in a diffusion problem (7). Allen and Severn (8) 
treated essentially the same assumption in a different way using relaxation 
techniques 

This paper describes two new ways of dealing with the general problem 
of a moving boundary. In the first method the boundary is rendered 
stationary by an appropriate change of variable. This means that instead 
of solving the usual diffusion or heat-conduction equation with a moving 
boundary we have essentially an eigenvalue problem with fixed boundaries. 
The equation to be solved contains parameters, associated with the move- 
ment of the boundary, for which values have to be determined such that 
the boundary conditions are satisfied. In the second method Lagrangian 
interpolation formulae for non-equal intervals are introduced into finite- 
difference formulae in order to follow the movement of the boundary. 


2. The mathematical problem 


The following particular problem will serve as an introduction. Suppose 


T 
( 





iat an infinite sheet of uniform material of thickness 2a is placed sym- 
metrically in a solution of extent 2/, and that the solute is allowed to diffuse 
into the sheet, which is initially free of solute. The sheet contains a number 
of fixed sites on each of which one diffusing molecule can be instantaneously 
and irreversibly immobilized and so withdrawn from the diffusion process. 
There are S sites per unit volume of the sheet. The concentration of solute 
in the solution is always uniform and is initially C,,, expressed as the number 
of molecules of solute per unit volume. The concentration at any time just 
within the surface of the sheet is taken to be that in the solution. If we 
denote by C the concentration of freely diffusing molecules at any point 
in the sheet in number of molecules per unit volume, the concentration of 
molecules immobilized on the sites is zero if C is zero, and equal to S when 
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C is positive and non-zero, because the immobilization is rapid and irrever-| while 
sible. 

If we take x to be the distance coordinate measured in the direction of} 
diffusion from x = 0 at one surface of the sheet and use the variables 


A = 2/a, T = Dt/a*, c= C/Cz, s= S/Cpz, (1) 


ee es At t 
where D is the diffusion coefficient assumed constant, we have to solve 
the equation ne 
oc O*C , r 
—=—, 0< X< Xj, (2) 5 
oT oX? 
subject to the initial condition, considering only half the sheet for reasons 
of symmetry, aa [oe ‘ 
apne: c=0, T=0, 0<X<1, (3) 
and the boundary conditions -— 
c= 0, Roy > @ (4)° This 
dX, (5) 2 é 
. (Dd) 
‘yy ‘ y > 
— ie f 
dd s\oX ee a tl 
xX) A so 
a , i 1" ; 
c=] 7 | (et s) dX, A = @, rSé (6) 
0 
- we , ae ; whe’ 
Here X, is the value of X at which c = 0, and for X <_X, all sites are | (16) 
occupied and for X > X, they are all unoccupied. Thus X, is a function 
of 7 and increases from zero as the boundary between occupied and free 
sites moves into the sheet. 
Ag 
3. Transformations (9, 
There are two difficulties in the way of treating this problem by finite- for ; 
difference methods. The first is due to the singularity on the surface of the i 
sheet, X = 0, and the second to the moving boundary, the position of 
which at any time has to be determined as part of the solution. The singu- whe 


larity can be handled by introducing new variables given by 
¢= 2/T", = 7, (7 


when (2), (3), (4), (5) become 


pel , — Bonk (8) 
0& Or c€ 
c = 0, & OC - 0 (9) 
c= 0, fF > €. 7 > 0, (10) 


where & 





Tevel 


10n oj 


les 


‘ASONS 


inite 


yf the 


no 





MOVING-BOUNDARY 





PROBLEMS IN DIFFUSION AND HEAT FLOW 223 


while (6) becomes 


c= 1—<- | (c+s)rdé é=0, 7>0 (12) 
dc de 
+ = 0 we have Gee aici 13 
At we have rT Te (15) 
( I é=0 (14) 
( 0 é €,, (15) 
(de > 
\aé}., 28h). (16) 


In deriving (16) from (11) we have assumed ¢€,/¢7 to be finite at 7 = 0. 
This is justified because we can assume the concentration on the surface, 
0, to remain constant for a very small time near t = 0, in which case 


v1 


the formal solution (9, p. 103) shows that X, varies as T", i.e. €, is constant. 


\ solution of (13) subject to (14) is 
c 1+-Berf(t€), (17) 


where B is a constant which together with &, is determined by (15) and 


(16) since we have | Berf(4é,) 0. (18) 
m£, erf(4é, exp($é,)? = 2/s. (19) 


\ graph showing the function on the left-hand side of (19) is available 
9, p. 105) from which a first estimate of the value of €, can be obtained 
for a given value of s. 


We now fix the moving boundary by writing 





£&), i 7 (20) 
when we obtain 
9 fC » OC  2/0c\ oc m 
oo a Se , 0< <i], (21) 
cy” CT S\CN} 1 7) 
| dé, | (| 4 €, n | (22) 
2 dr STE, cy 1 2r 
1 
a 
ae (c+s)r€, dy, 1 = 9, (23) 
0 


where we have used (¢c @7), to mean (éc/@n)¢_¢,. The problem is now that 


of solving (21) for trial values of €, and (éc/éy), which are mutually con- 
sistent with (22), until we obtain the solution for which (23) is satisfied. 








224 J. CRANK 


4. Finite differences 


There are doubtless several ways of obtaining solutions of (21), (22), and 


(23). We write 


—9/Ap 
iw -(<). (24 
s O71} 4 


and denote by c,, the value of c at the point m4» at time 7 and by c;* the 
corresponding value at time 7+-5r. Similarly, if A+, €°, and (éc/én);* are 
the values of A, €, and (éc/é@n), at r-+57, and nr is written for 7, a con- 
venient and workable method is based on the following finite-difference 
forms of the above equations: 


aT ;> , 1(2n+-1)f€2 + (€2)+(e+ —e 
\(8n) 407 | . 7 2( fi + (87) s( m al 
2 
x (dn)? (Cm +1 +Cm+1)— 2(Cm Cn) + Cu—2} T 
' (A 1A *)y 4: . 
zg 48 {cx iten 1 -Cm—1f> (25) 
ce 30 
cy 1 én ' 


We calculate values of c at equal intervals 5y in the range 0 < y < | at 
7 = 0 from (17), (18), (19), (20), and successive steps in the solution for 
the interval n 67 to (n+-1)67 are then as follows: 





(1) Estimate €/ at r = (n+1)6r7 and calculate (éc/é7);° from (26). 
(2) If we take m = M when » = 1, and use the approximation 
an \+ + 4 
c¢\" _ ¢-2—4¢m-1 (21) 
On), 20 
since ¢,, = 0, with the special case of (25) for which m M—1, we have 


two equations from which to calculate c,,_, and ¢,,_». 

(3) Integrate backwards along the line tr = (n+-1) 87, using (25) to obtain 
c;,-, knowing c;, and c;,, for successive values of m. 

(4) Adjust the estimate of £} until the value of c at 7 = 0 so obtained 
agrees with that calculated from (23), where the integral is obtained by any 
convenient formula for numerical integration. Continue in successive steps 
dr until X, = 7€, = 1 when this part of the solution ceases since all the 
sites throughout the sheet are then occupied. To continue further we merely 
have to solve the simple diffusion equation with no complications associated 
with a moving boundary. 


That this is a powerful method of dealing with this relatively difficult 








mMovt. 


proble 


and s€ 


Sint 
and t] 
theref 
when 
concel 
centre 


of col 


oi mo 
imply 
out. | 
tion 0 
temp 


e.g. W 


consi 
fluid 
into t 
or lib 
as an 
be ay 
If we 
the n 


5092 





MOVING-BOUNDARY PROBLEMS IN DIFFUSION AND HEAT FLOW 225 





problem is illustrated by the results obtained for the case l/a 1,s = 0-25 


-and and set out in Table 1 below. 


TABLE | 








24 
Lia l. Ss 0-25 
; \ i €, are tabulated at r = 0-5 
Tie 
+ ar | fs f 
CON: ; 21 245 17 108 50 ° 1*802 
rence 
244 108 °) 1812 
4¢ 108 I 
Since X &,7, we see that when 7 = 0-5 and €, 1-812, X, = 0-906 
nd the boundary is 3ths the way to the centre of the sheet. Table 1 shows 
therefore that in this case quite reasonable accuracy is obtainable even 
en only two steps in time and four in distance are taken to evaluate the 
(25 meentration distribution at the time when the boundary has reached the 
ntre of the sheet. The number of steps in 7 and 7 depends to some extent, 
PY f course, on the parameters //a and s, and on the accuracy required. 
5. Amore general problem 
1 at The initial condition which we have considered is the one likely to be 
mn. Tor ost interest in diffusion problems. In a heat-flow problem it would 
ply that the sheet is initially at the transformation temperature through- 
But a more general problem is likely to arise in which the transforma- 
m occurs at some temperature intermediate between the initial uniform 
mperature of the sheet and that to which its surface is raised at zero time, 
27 g. we may have a sheet of ice initially at —10° C. whose surface is raised 
(, In such a case we have to consider the flow of heat on both sides 
have fthe boundary surface on which the transformation is occurring at any 
remembering that the heat diffusivity may be different on the two 
btair S 
The problem in heat-flow corresponding to the diffusion problem we have 
ained msidered in previous sections is that of a sheet immersed in a well-stirred 
y an} id or a perfect conductor of limited extent. As heat flows from the fluid 
steps to the sheet a transformation boundary, on which latent heat is absorbed 
Hi the t liberated, moves into the sheet. We shall continue to use this problem 
1erel} is an illustration, though it will be clear that the method of solution can 
ejated be applied to any known boundary condition at the surface of the sheet. 


Tf ye ° . ~ 
ii we denote by suffix J properties between the surface of the sheet and 


fficult the moving boundary, and by suffix JJ those between the boundary and 
, Q 
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the centre of the sheet, the relevant equations in terms of the variables , 


and 7;, where 7, = Tt = (D,t/a*)t = {K,t/h, a, (25) 
76 sie Oe 2y 06,( (06, 06.,; 
are $——! ox f87,—! 4 : u( ) - K( = \ Oe ye, (29 
ep ér,* 8s en\\enh en hy) 
076 »o of 2 06 fale) (og 
2p = fie, 4 Z: | oe an oé | ge (30 
cy” CT] s on | C7) 1 cy J 
l dé, ; l (5!) ; K(f) - E, Bn) 
2dr, 87, ,| on), on 1] 277 
1 r 
h aft : a fh : 
°(1—8,) (0,+8)7,&, dy | 16,,7,&,dy, 1n=0, (32 
h, . ty ky, 
0 
6, 61 6, n l 33 
O77 Q, 7) 00, (34 


provided the sheet is still effectively semi-infinite, i.e. the temperature at 
its centre has not changed appreciably from its initial value to within the 
accuracy of working. In the above equations @ is temperature expressed 


as a fraction of @,, the initial temperature of the well-stirred fluid, and | 


measured on a scale such that the sheet is initially at zero temperatur 


throughout; A = K,,/A,, the ratio of the heat conductivities on the tw 
sides of the moving boundary; D = D,,/D,, the ratio of the heat diffusivi- 


ties: hy, h;, hj; ave respectively the heat capacities per unit volume of the 
well-stirred fluid in which the sheet is immersed and of the two regions of 
the sheet. The relationship D, = K,/h, is implicit in (28). Also 

8 L/(o, 9p), 


where L is the latent heat absorbed per unit mass of substance undergoin 
the transformation and a; is the specific heat of substance J. The solutior 


at T, O1s 6, col B,erf(3é, ”)s () - n - l. (35 
Gi; 7 oe ee 9 
and 6,, = Op —- at ert £1 1 erp $1|, (36 
l—erf{é,/(2D')}| 2D! 2D} 
where &, is given by 
} ( wee 
(1—6,,) o-étla_ KO, e781 —s ar 
7 erf(4é,) D'7* erfef{é, /(2D*)} alli 
and B, by 1+ B,erf(4é,) = Op. (38 


Writing the equations in finite-difference forms as before and evaluating 
6, and 6,, at tr; = 0, we find that a convenient procedure is: 

1, Estimate €;° at +, = (n+-1)67, and calculate (€6,/6y) 7 — K(€0,, )i 
from the finite-difference equivalent of (31). 
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2, Estimate (€6,,/€), and find by solving (30) the value for which 0,, = 0 


when 7 XD. 

}. For this value of (€@,,/¢y), find the corresponding value of (€0,/éy), 
onsistent with the estimated €; , and integrate (29) backwards from 7 = 1 
to find the value of 07 at 0. 

t+, Take new estimates of €; till the value of 07 so obtained agrees with 
hat caleulated from (32). Continue in successive steps 7, until a change 

the temperature at the centre of the sheet occurs which is appreciable 


to within the accuracy of working. 


6. Two transformations 

When the temperature starts to change at the centre of the sheet the 
X l, 
€ &, has been used to fix the boundary at X 


symmetry condition ¢6,,,cX = 0, must be taken into account. 


X/X, x: 


which the transformation occurs, but the central plane of the sheet, 


The variable 7 


\ 
gested in conversation that it might be practicable to fix both the boundary 
ik=a2 
region I] 


suggest themselves 


1, does not occur at a constant value of 7. Professor Hartree sug- 
] YP 


and also the central plane by using a second space variable in 


[f we continue to work in the (€,7) plane the following variables 


c. Vv ‘9 Sas (39) 
y $1 < 
C= 2, & SE <1, (40) 
1—§y 
e € t, at the centre of the sheet, X 1. The relevant equations 
ecome 
-H cou ob dé 
I ¢2 T{ ¢2 >. wae 
- S17] (3 f1719 \, 0 " I (41) 
; atid a OT] 
728A cA R . dé cA n 
Dr je Seen a ry E:)(E "14,) w @<ct<t 
ee ( aT] C4 
(42) 
¢ ! ¢ 
| dé, 1 {] (©) Kr, ea \ & (43) 
2 dr; 77 8\Sq\E J 1—£,7,\ €C J) 2T) 
; 0,=0,= 07, y=! (=60 (44) 
06,,/0€C (), C ] (45) 
1 1 
h a F a i h 
0 ai - v . 
(0,+-8)7,&, dr (l—7,£,)0,,df, 7 = 90. (46) 
/ Pore J hy , 
0 0 


lhese are of the same form as (29)-(34) and can be treated in the same way, 


the computational labour involved being only slightly heavier. 
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7. Lagrange interpolation formulae 

An alternative is to return to the (X, 7’) plane and to extend the finite. 
difference formulae which are commonly used as approximations to deriva. | 
tives in such a way that they can be used to follow the motion of the 
boundary on which the transformation occurs. This means developing 
finite-difference approximations to derivatives based on functional values 
which are not necessarily equally spaced in the argument. If we make use 
of the Lagrangian interpolation formula 








f(x) > I(x) f(a;), (47 
{=0 
where l(a) — PrlX) (48 
: (x A ;)Pp(;) 
and Pp, (2) (v7—A,) (x a,)...(x -An—1)(%—Ay), (49 
and restrict ourselves to three-point formulae, i.e. n = 2, then we find 
1 d*f(x) f(a) f(a) f(a) - 
-- = Xl —_—— + —— ———, () 
2 dx (4g —4,)(Mg— 4g) (4, —Ag)(4,—Ag) — (€g—Ag) (4g—) 
d (: ° . » . , » 
and je = 19(2) f(a) +1 (x) f (ay) +-13(2) f (a2), (51) 3 
ar 
where 
ia ot SC. 
(4) —4,)(Ayp—Av) (a, —@,)(a, —y) 


Bleed an (52 
(@,—Ap)(4g—A}) 

We apply these formulae in the neighbourhood of the moving boundary 
Consider the sheet to be comprised of .V layers each of thickness 6X and 
let the boundary at time 7, be somewhere in the (7+ 1)th layer such that 
X, = (r—1+p)6X, where p is fractional and 1 < p < 2 as in Fig. | 
This is found in practice to be preferable to the apparently simpler ex- 
pression X, = (r+p)dX because of difficulties associated with p <1. 
Then if we take, for X < X,, 


F (4) 0-2 f(a) O,1 f (42) On, (53 


the expression (50) becomes 


( "6, f 2 0, 2 0, 1 Oy 


: I 
| 1)} 


ox?) ene (SX)2|p -i p Tl p 
and from (51) we have 


( 3) - ] { p@,. 2 Pa a) A. 2p +1 a,,.\ (55 
X=X) 


~ 8X |p+l "ppl 77 
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? r+] r+Z 

{ Qin Qa, a> 

Fic. | 
If similarly we take for X As, 
f(a,) 6 f(a,) 6.45, JlG@e) = 6,12; (56) 


we find 


2 6, A. . 
4 = — — j frie) (57) 
OX? } y_w@ansx dX )?|(p—2)\(p—3) 2—p' 3—p| 
0 2n 5 3 2 
(7) A Ee a EE 6, .0\. (58) 
X \ dX |(p—2)(p—3) p—2 3I—p “| 


On combining these formulae with approximations of the Schmidt type, i.e. 


0: H oH ? a ~ 
f me. — ot tame 7’, (599) 
57 C 
we have finally 
207,10... 8, Op 
i <0 olf Et... ae | (60) 
(OX )*|p+1 P p(p+1)J 
257’, bn 6. Ca 
| 0. 1} 1 r+i r+2 | (61) 


5X *\(p 2)(p a) £ , 3 pl 
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These are to be used in conjunction with the usual Schmidt finite-difference 
formulae 


2 é-+ ax (6,.41—26,,+6,-1) lamgr—2, (63 
6- = 0.4 Sry (Ome 26,,+9m1)) T+2<m< N—-1, (64 
0% = Oy ses i—9y), X=], (65 
(05 A4)( ba 2lho 7 : 267; (0,—9), X 0. (66 

ah, 5X (6X)? 


The surface condition (66) comes by imagining the sheet to be extended 
one step 5X beyond the surface. The equation of the type (63) for 4, is 
then combined with a finite-difference form of 

fal - 06 

h,!- K,—, x= 0, (67 

ct Cx 
expressing conservation of heat, in order to eliminate the hypothetical 
quantity @_, (9, p. 198). When the boundary surface passes from one 
layer to the next, i.e. 7 increases by unity, we obtain the value of @, by 
using the Lagrangian interpolation formula 
l—» 2(p—1 2 

Pp (2 fa Ds 


—27T r-1T 


_ @ 
lt+p ' p p(p+1) 


T: (68 
and this becomes the new 6,_, to be used in (60) and (62). 

The essential difference in using the two methods described is that the 
application of Lagrangian interpolation formulae involves many steps in 


time but each step is relatively simple to work and the computer does not ! 
have to exercise any judgement. When the boundaries are fixed by using . 


suitable variables, far fewer steps in time are needed but each may involve 
several trial solutions. The number of trials needed depends to a great 
extent on the ability of the computer to make good estimates. 
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SUMMARY 
Initial-value systems, particularly involving first-order differential equations, can 
be transformed into systems of higher order, and treated either as initial-value 
problems or as boundary-value problems. Two such transformations are considered 
in detail, and the advantages and disadvantages of these methods are discussed and 
illustrated with particular reference to linear systems. 


Introduction 

1. An ordinary differential equation with all relevant boundary condi- 
tions prescribed at one point gives rise to an initial-value problem, and its 
solution is usually effected by a step-by-step method. Boundary conditions 
prescribed at two or more points result in a boundary-value problem, and 
different techniques are commonly used for its solution. 

Allen and Severn (1) suggested a method whereby a first-order system, 
normally solved by initial-value techniques, might be replaced by a second- 


order boundary-value problem. No indication was given of the merit of 


this procedure, except when extended to certain partial differential 
equations. Fox (2) also proposed to treat a derived second-order equation, 
but remarked that either initial-value or boundary-value methods could 
be used, the choice depending on the nature of the solution. No great 
details were given regarding this choice, the paper being concerned mainly 
with the selection of the most accurate simple finite-difference equations. 

The purpose of this paper is to compare and examine these methods in 
more detail, and to indicate their fields of useful application for linear 
differential equations. 

We show that boundary-value methods have advantages only when 
initial-value techniques are unstable, in the sense that rounding errors 
accumulate catastrophically and swamp the required solution, and that 
Fox’s method is then likely to be superior to that of Allen and Severn. 
Boundary-value techniques can sometimes be used, but with caution and 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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discretion, when step-by-step methods are stable. Caution is needed since 
the differentiation of the original problem may introduce extra solutions 
whose effect is to make the boundary-value methods lack precision in virtue 
f the near-singularity of the matrix of the finite-difference equations. 
The introduction of extra solutions may also affect adversely the use of 
nitial-value techniques for the derived problem of higher order, and these 
methods are successful only when the unwanted solutions decrease more 
rapidly, or increase less rapidly, than any essential component of the 


solution of the original equation. 


The two methods. General comparison 


2. For the solution of a linear first-order system of the form 
y +f(a)y g(x), y Yo at x ce (1) 


\llen and Severn (1) introduce an auxiliary function z, connected with y 
vy the equation , ‘ 
‘ 1 y k(z’,z). (2) 
The form of k(z’,z) is to some extent arbitrary, but is generally chosen to 
be a simple linear combination of z’ and z such that substitution in (1) pro- 
luces a second-order equation for z in which the coefficient of z’ vanishes. 
For example, if ‘a , 

y = 2'—f(x)z, (3) 
substitution in (1) gives 


(f'+f?)z = 9. (4) 


le suggested boundary conditions to be used with (4) are given by 


fz Yos z Z, (5) 


<0 o°,U n 
where the first condition of (5) is the equivalent of the initial condition 
for the original first-order equation, and Z is an arbitrary value of z at any 
other point 2, in the range. It is easy to see that this value can be arbitrary 
without affecting the resulting solution for y, which is recovered from 
3) when (4) and (5) have been solved for z. 

\llen and Severn suggested that 2, should be a point in the range of x 
beyond which the solution is not required. It is, of course, equally possible 
to choose x, to be 2, giving an initial-value second-order problem. 

3. Fox (2) produces a second-order differential equation by direct 
differentiation of the original first-order equation and subsequent elimina- 


tion of the first derivative. For example, (1) is replaced by 
y"+(f'—f? ly = 9 —S9. (6) 


For boundary conditions he takes the given condition y = y, at x = 2p, 


ind the first-order equation (1) applied at 2, or x,, according to the 
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particular technique, of initial-value or boundary-value type, which }; 
wishes to use. 

4. Mathematically there is little difference between the two methods | 
Numerically, however, Fox’s method is superior in the sense that it permits| 
the use of simple finite-difference equations of very small truncation erro; 
Both methods have this facility as far as concerns their respective second-| 
order differential equations, for which three-term formulae can be estab.| 
lished whose major truncation error is 6°y/240. For the boundary condition) 
which involves a derivative, however, Fox can use a three-term formul 
whose main error term is .5°y/90, while the corresponding equation for th 
other method, the first of (5), has in its finite-difference form of equa 
simplicity an error as large as .5*y/3. The reason for this is connected wit! 
the fact that Fox’s boundary condition (1) of this type holds not only at th 
boundary but at all other points also, whereas the first condition of (5 
holds only at the single point 25. 

Since the first-order system can be solved by using a two-term recurrene 
formula whose main error is 6°y/12, the initial-value method of Allen and 
Severn can never be advantageous. Fox’s initial-value method would 
appear to have no advantage over a process which solves the first-order } 
system with a three-term formula whose main error is d°y/90, but this 
process suffers from the occurrence of unpleasant oscillation in the pivota 
values, as illustrated by Fox (2). Both the new initial-value methods 
remove this oscillation, but only that of Fox preserves both accuracy and 
simplicity in the finite-difference formulae. 


5. A small truncation error, however, is not a sufficient criterion on 


which to base a claim for a particular method. Stability is also necessary 
and this we proceed to discuss. In particular, the lack of stability in 
initial-value methods must be the only reason for using boundary-value 
techniques, when there is a choice, since the former are in principle easier 


to apply, and the truncation errors are comparable in both techniques. 


Stability of boundary-value techniques } 


6. All boundary-value techniques involve the solution of a set of simul- 
taneous algebraic equations, linear if the differential system is linear. If 
the equations are well-conditioned it is easy to get a solution, either by 
direct methods or by iteration or relaxation, and small changes in the 
coefficients have little effect on the solution. With ill-conditioned equations 
the associated matrix is almost singular, and small residual errors may 
have a large effect on the solution. If the matrix is truly singular it means 
that the homogeneous equations, obtained by replacing the constant terms 


by zero, have a non-trivial solution, an ‘eigensolution’, and the original , 





problen 
constan 
This 
differen 
In part 
eigenva 
probler 
initial - 
ifour d 
or neal 
solutio: 
legitim 
eS 
and th 
the de 
Con 


of whi 


given 


The d 


and t] 


Subst 
of A | 
If x, 
(10) i 
Aist 
be sa’ 
becor 
8. 
well 1 
has i 
to th 
the f 








ich | 





BOUNDARY-VALUE TECHNIQUES 235 


mroblem either has no solution at all or an infinity of solutions, one arbitrary 
eonstant being undetermined 

This situation, of course, may arise quite naturally from the finite- 
lifference equations used to represent an arbitrary second-order system. 
In particular the determination of such a solution, or of a parameter, the 
eigenvalue, which permits such a solution, is a respectable and well-defined 
problem, and no question of stability arises. For a differential system of 
nitial-value type, however, there is no possibility of an eigensolution, and 
four derived second-order boundary-value problem gives rise to singularity 
or near-singularity of the finite-difference matrix the cause is the extra 
solution introduced by the differentiation. The epithet ‘unstable’ can then 


egitimately be attached to the method. 


7. It is easy to show, by simple examples, that both the method of Fox 
ind that of Allen and Severn will often be unstable in this sense, though 


the degree of instability is controllable to some extent. 


Consider, for example, the differential system 


y jt 0, y(O) l, (7) 

f which the solution is y |. Fox’s derived boundary-value system is 
given by 

y y | 0 y(0) a v(x.) y(x,,) ] 0. (8) 


The differential equation in (8) has the general solution 


y Aet+ Be-*-+-1, (9) 
ind the use of the first boundary condition gives B —A, so that 
/ A(e*—e-7)+-1. (10) 


Substitution in the second boundary condition permits the determination 


tA from the equation > Ae—*n 0. (11) 


If x, is «©, (11) is satisfied for any finite value of A, and Fox’s solution 
l0) is then different from that of the original problem. For any finite ~,, 
1 is theoretically zero, but from a numerical point of view equation (11) can 
be satisfied approximately with finite values of A. The discrepancy in (11) 
becomes smaller, for fixed A, as x, increases. 

8. At an interval for which the complementary functions are reasonably 
well represented by the chosen finite-difference equations, this phenomenon 
has its counterpart in the difficulty of obtaining precise numerical solutions 
to the simultaneous finite-difference equations. For example, if we replace 


the first of equations (8) by the simple general equation 


Yr4y—(2+h?)y,+-y,. +h? 0, (12) 
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the third equation of (8) by 
, Yn+1 -Yn-1 2hy,+ 2h ‘ Q, (13 
and eliminate y,,,, from (13) and the particular equation of the set (12) for 


which r = n, we have, with interval h = 0-2 and n = 10so that x, = 2. 


), 
the simultaneous equations 
y(0-0)—1-00 = 0 


y(0-0)— 2-04y(0-2) f y(0-4) +()-04 — 0 
y(0-2)—2-04y(0-4)+- y(0-6)+ 0-04 0 


y(1-6)—2-04y(1-8)+-y(2-0)+ 0-04 0 





y(1-8)—0-82y(2-0)—0-18 0 


From these equations, solved for example by the relaxation method, it is , 


possible to obtain a two-decimal solution, shown in Table 1, which is far 
from the truth though the residuals are apparently negligible. 


TABLE | 


oo o°2 o*4 06 08 I°o 1*2 I*4 1°6 18 2"0 
I°co T'O4 1°08 I°l2 118 1°23 1°30 1°38 1°48 I°5Q I*72 
K 0°00 0°00 oOo! ool foe ove} | 0"00 0*00 oO'ol 0°00 0°00 


The correct solution of (14) is that of the differential equation, y = | 
everywhere, and there is no great difficulty in obtaining this solution. In 
many practical problems, however, the coefficients may have to be rounded, 
and the rounding may cause significant errors. Moreover, many guarding 
figures would be needed in a direct method of solution, while relaxation 
or iteration would be slow and uncertain. This is a definite example of 
instability. 

Fox (2, p. 376) noted a similar phenomenon in a different example, but 
the ill-conditioning was less severe because the total range of x was smaller. 
This confirms an early comment in section 7 and suggests that limitation 
of the total range, if possible, may reduce the degree of instability. 

9. For the example of section 8 the finite-difference equations are of 
simple form, with constant coefficients, and we can obtain, by known 
methods, an analytical solution which is interesting and which is needed 
in later investigations. To illustrate instability we assume that the main 
finite-difference equations (12) are satisfied exactly, but that the second 
boundary condition, the third of (8), is not quite satisfied, and we represent 
this by adding a term 2e to the right of (13). We then find the solution 


esinh ré 
sinh 6 cosh n@—h sinh n@’ 


where cosh 6 = 1+ $h?. (16) 


Y, (15) 
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From (16) we see that sinh @ = h, if h is small, so that for any appreciable 
yalue of n the denominator in (15) is small, and the error term may be 


substantial even for small «. With h = 0-2 and n = 10, for example, 
sinh @—/, 0-0010, and the error in (15) is about 30esinhré. 
The smallness of the denominator corresponds, of course, to near- 


singularity of the finite-difference matrix and its effect, as before, is reduced 
DY restricting the size of n. 

10. The method of Allen and Severn also has instability defects. For the 

ex umple of section 7 they use the equations 
0 2'(0)+-2(0) 1, Z Z, (17) 

ecover y from the equ ition 

y 2’ +2, (18) 
lhe simplest finite-difference form of the first of (17) is identical with that 


of (12), with y of that equation replaced by z, and the second equation 


2h. (19) 


f (17) is replaced by 25 —2_1 + 2hezy 


If we add an ‘error’ of amount 2e to the right of (19), we find for z, the 


solution 
? sinh 6 cosh 76 —A sinh ré e sinh(n r)@ 
1+(Z—1) ! on} =p-pasoraciemione, f00) 
sinh 6 cosh n@—A sinh né © sinh @ cosh né—A sinh nd 
which is numerically imprecise for the reasons given in section 9. The 
luation of y. from 
2) r ~ )},- » 
2hy 21 — 2p—-1 + 2hz,, (21) 
the simplest finite-difference replacement of (18), is obtained from 
Z—1)h* sinh rd e sinh 6 cosh(n—r)@—A sinh(n—r)é 
t(sinh @ cosh né sinh né) = h sinh 6 cosh n#—h sinh né 
(22) 
{which at least the last term represents a substantial error, again greatest 
r larges 
11. We note that replacing e by zero in Fox’s result (15) gives the correct 
solution of the differential equation, whereas the method of Allen and 
Severn has an error of amount depending on Z. This quantity does not 
ippear, of course, in the solution of the differential system, and its presence 


in the finite-difference result for y is a consequence of the use of approximate 
inite-difference equations. This point is also important in later investi- 


rations. 


Stability of initial-value techniques 
12. We turn now to an examination of initial-value techniques from the 
point of view of stability. The first-order system (1) has the solution 


y AOQ(x)+ P(x), 
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in which C(x) is the complementary function, the solution of the homo. 


geneous form of (1), P(a#) is a particular integral, and A is a constant 
determined from the initial condition. 

The best numerical method of solving the first-order system directly 
uses a two-term finite-difference formula in a step-by-step manner, and 
the solution can be written in the form 


y, = AC,+P,, (24 
analogous to that of (23) for the differential system. Apart from truncation 
errors, rounding errors are propagated in amounts proportional to the 
‘complementary function’ C,. If C, decreases faster than P. these errors 
have negligible effect, and the process is stable and gives good results. If 
C, increases faster than P, the error also increases and, if the possible 
rounding error at any point x, is e,, the solution obtained is incorrect by 
an amount of the order of «,C,/C,, where r > s. If A is large this error is 
unlikely to be serious, since the number of correct significant figures 
obtained is approximately constant, and this is the likely requirement. 
The error is serious, however, if A is small and in particular if A = 0, 
For then the solution, which should consist solely of the particular integral, 
is ultimately swamped by the accumulation of rounding error. 

Consider, for example, the solution of the system (7) with the assistance 
of the finite-difference equation 


Yriy( l—th)—y,(1+$h) h. (25 
The general solution of (25) is 
2+h\" 
y= 1-4 Al ‘ (26 
2—h 
The initial condition gives A 0, but if we commit an error by taking 
Yo = 1+e, a circumstance differing little from the commission of a rounding 
error at the next step, we find A e and 
| h r 
- ‘ie 
Y, 1 (5 i} . (27) 


of which the second term is likely to increase rapidly for ‘reasonable’ 
values of A. 

13. Initial-value techniques used with the derived second-order systems 
offer no advantage in this case. For example Fox’s derived problem, with 
the simplest finite-difference formula and with an error introduced into 
the initial condition which involves the first derivative, has the solution 

y, = 1+. cosech 6 sinh r, (28 


and the second term grows without limit. 


Indeed, the introduction of a second solution, through differentiation, 
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ym may introduce instability even when the corresponding first-order process 
tant s stable. For example, the system 

4 y'+y—1 0, y(0) l (29) 

- (2—h\" 
2s ne Ss j t t = 9 30 
an has the table solution /, l (5 i) (3 ) 


ysing the method and notation of section 12. Fox’s derived initial-value 


system 1S given by 


+) ? 

. y”"—y+1 = 0, y(0) = 1, y'(0)+y(0)—1 = 0, (31) 

~~ ind if we introduce errors « and 27 in the value of y(0) and in the finite- 
i lifference form of the last of (31) respectively, we find the solution 

sibl y, | + e cosech 6(sinh 0 cosh ré—h sinh ré)+ 7 cosech @sinhré, (32) 

CDY y of which the last two terms increase steadily with r. 

ee It would appear that solution of Fox’s derived second-order initial-value 

rure 


system has advantage over that of the original system in one case only, 
when the complementary function of the latter contributes in substantial 
measure to the true solution and when the differentiation has introduced 


Br in extra solution of smaller magnitude. This would be true, for example, 


for the differential equation of (7) whose boundary condition, say y(0) 0, 
- illows the complementary function e* to appear in the solution: the 
differentiation here produces the ‘less important’ function e~’. 
Boundary-value techniques for unstable initial-value problems 
os 14. So far we have suggested no method for the effective and accurate 
‘ solution of problems like (7), for which any step-by-step method is unstable. 
king We find, however, and proceed to demonstrate, that the boundary-value 
ding method of Fox and, less conveniently, that of Allen and Severn, can be 
used efficiently for this purpose. 
We recall that the boundary-value techniques, for solution of the second- 
21 order problem derived from (7), are defective only through the introduction 
,;,.° of near eigensolutions of the corresponding homogeneous problem. The 
oe finite-difference equations are likely to be affected in an analogous manner 
only if the chosen interval / is small enough to permit adequate representa- 
Ue tion of the complementary functions. If we take an interval large enough 
wi to exclude adequate representation of these unwanted functions, but small 
y enough for good representation of the required particular integral, we are 
101 likely to avoid the boundary-value type of instability and obtain a good 
28 solution 


his conjecture is verified in the case of Fox’s boundary-value solution 
15) for the problem of equations (7). For if A is so large that sinh 6 differs 
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appreciably from A the denominator in (15) is no longer very small, we cay 
approximate to sinh n@ and cosh n@, for reasonably large n, by e”?, and the 
error term in (15) is everywhere less than ¢/(sinh 6—h). 

We conclude that. Fox’s boundary-value method will give good results, 
when step-by-step solutions are unstable, if the finite-difference interval 
is chosen with sufficient care. 

15. The boundary-value method of Allen and Severn, however, has q 
severe defect. Although the condition on y may exclude the complementary 
function from the solution for y, the auxiliary dependent function z contains 
this complementary function if the quantity Z is arbitrary. The numerical 
differentiation involved in the recovery of y does not then completely 
exclude Z, as we observed in the finite-difference solution (22) of the method 
of Allen and Severn for the problem of equation (7). Ifh is so large that the 
denominator in (22) is not small, the term in e is effectively negligible, as 
in Fox’s method. There remains, however, a truncation error of approxi- 
mate amount (Z he 


f= — 
4(sinh 6— hyen—6 "0 


(33) 
and for large h this may be quite substantial. 

There are two evident ways of reducing this error. The first is to choose 
Z, in this case equal to unity, so that the large complementary function is 
excluded from z,, and therefore from y,. 


r r 


In more complicated problems 
it is difficult to see how this can be achieved without considerable numerical 
labour. The second method is to choose the value of n, and hence the total 
range, so large that for the values of r needed the denominator in (33) be- 
comes large and hence 7, small. This is more practicable, but again involves 
more work than is needed in Fox’s method. 

16. It might be thought that step-by-step methods should also giv 
good results if the interval is large, so that the complementary functions 
are inadequately represented. Examination of our finite-difference solu- 
tions, however, reveals the limitations of this suggestion. 

For the step-by-step solution of the derived second-order problem, for 
example, we see that the error term in (28) increases steadily with h, and 
this method is always unstable. 

We can apparently do better with the first-order method, since increasing 
h, beyond h 2, ultimately decreases the error term in (27). A disadvan- 
tage of such a large h is the fact that the complementary function, ceasing 
to have any connexion with that of the differential equation, here becomes 
an oscillatory function, whose adverse effects have already been noted by | 
Fox (2) in other contexts. With such a large h, moreover, the required | 
particular integral may cease to be well represented in wensnileetl problems 
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Fox's boundary-value method has an apparent advantage for ‘reasonable’ 


}. With h = 1, for example, the error term in (15) is approximately 
gee"? nowhere exceeding 8e, and only one ‘guarding’ figure is needed: 


with this h the error term in (27) is 3’e, which grows rapidly. 

17. Problems in which the complementary functions are completely 
excluded by the initial conditions may be thought artificial, and of some- 
what rare occurrence. They arise quite naturally, however, in the tabulation 
of certain mathematical functions for large values of the argument. 

For example, Fox and Miller (3) showed that the exponential integral 
could be effectively tabulated for large 2 by using the auxiliary function 


T(z) ez Ki(1/z)—z, (34) 

where 7’ satisfies the first-order system 
T’ —z?T +1 0, T(0) = 0, (35) 
differentiation being with respect to z = 2-1. The solution required is in 


fact the particular integral, recognizable in the form of the asymptotic 


series ‘i g2@+ 212843! 244+... (36) 


The complementary function of the first-order equation is e~'”, and in the 
more important negative range, say from z = 0 to —0-1, corresponding to 
rranging from 00 to 10, direct step-by-step solution of (35) is impractic- 
ible, either through catastrophic accumulation of rounding errors, with 
small h, or through the type of oscillation mentioned in section 16, which 
occurs for relatively large A of the order of 0-005 or 0-01. 
We can, however, derive Fox’s second-order system 
I" —(z+—2z)T+2T = 0, T(0) = 1, T’ (2,)—z,? T(z,) +1 = 0 
(57) 


ind although step-by-step methods, corresponding to z 0, are quite 


n 


unstable, the boundary-value technique, used say with z, = 0-1 and 


h= 0-01, is quite stable and easy to apply. Fox and Miller did not use 
the third of (37) as a second boundary condition, but took a value of 
['\0-1) obtained from previous calculations. The use of the third of 

17) does not affect adversely the condition of the linear equations (it is, 

lact, somewhat improved), and dispenses with the need for other 
computational assistance. 

Though we have concentrated on the solution of first-order systems, 
similar methods might be used for initial-value problems of higher orders 
ior which step-by-step methods are unstable. For example, the system 

y” — lly’ —12y+ 22e? 0, iO) = I, y'(0) ] (38) 
has the solution y = e*, but the complementary solution contains e!*”, and 


“2 R 
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rounding errors increase catastrophically in any step-by-step method. The 
problem can be solved, however, and in a variety of ways, by transforming 
it to boundary-value form. 


Conclusion 


18. We have investigated two methods by which first-order differentia 


systems can be transformed into second-order systems, soluble either by 
initial-value or boundary-value techniques. The illustrative examples 
used to indicate the fields of application of these various methods, have al 
involved linear systems with constant coefficients, for which the finite. 
difference equations have known analytical solutions. For these problem:' 
we have reached the following conclusions. 

When the two-term step-by-step solution of the first-order system i 
stable: 

(i) the second-order system of Allen and Severn, solved by step-by-stey 
methods, has no advantage ; 

(ii) the step-by-step solution of Fox’s second-order system may give 
more accurate results, for the same labour, in virtue of its more accurat 
finite-difference equations, and the fact that they do not give rise t 
unpleasant oscillation (we have not used the most accurate finite-differenc: 
equations in the illustrative examples, since they do not affect the questio 
of stability). This is true, however, only when the differentiation involve 
in producing the second-order system does not introduce a solution whic! 
grows faster than any solution of the original problem ; 

(iii) boundary-value methods applied to the derived system may b 
unstable through the near-existence of eigensolutions. When this can be 
controlled by limiting the size of the total range, the smaller truncatior 
error in Fox’s derived system may give it some advantage over all other 
methods. 

When all step-by-step methods applied to the original system ar 
unstable: 

(iv) boundary-value techniques applied to the derived systems may giv’ 
good results, and in this case the method of Fox is more convenient tha! 
that of Allen and Severn ; 

(v) similar considerations also apply to unstable initial-value systems 
of orders greater than the first. 

Some of these conclusions, of course, may need modification in detail fo 
more complicated problems. We note, for example, that the two derived 
homogeneous systems are not identical when the original equation ha: 


» . *y: y . at 
non-constant coefficients, and this may affect stability. Non-linear equa 


tions, also, will involve other considerations. We have, however, indicateé 
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the reasons for the conclusions reached, and similar considerations may 
guide a computer to the best choice of method in a given practical problem. 

Finally, though our conclusions are unfavourable to the method of Allen 
ind Severn, they in no wise criticize the use of this method for the solution 
f certain partial differential equations, the problem for which it was 
devised. 

We are grateful to the Director of the National Physical Laboratory for 
permission to publish this paper. 
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A FAMILY OF INTEGRALS OCCURRING IN THE 
THEORY OF WATER WAVES 


By E. M. WILSON (Admiralty Research Laboratory, Teddington) 


[Received 15 June 1956] 


SUMMARY 
An analysis is made of integrals of the form 


. ‘ cos 
| e2*8ec78 s6en® (Bsec 6) dé 
sin ° 


0 
directed towards finding efficient methods of computation. Particularly useful results 
are differential equations in the 8-direction and expansions in terms of Bessel functions 
which approximate to the integrals asymptotically. 


1. Introduction 


INTEGRALS of the form 
: . " COs : 
| e~*secV sec"? . (Bsec@) dé, (I 
; sin 
0 


where 7 is an integer, arise in discussions of various water-wave problems 


(cf. 1, 2, 3). The present paper describes an investigation of the integrals 
" : , 


paying attention particularly to properties which provide economic means 
of computation. We write 


le 


E" (x, B) E” (x, B)+7E" (a, B) M. etn — | e-esec?O sagng eibsecd dp. 


n 


0 
(2 
Functions of orders n 1, 0, 1, 2, and 3 are the ones of most immediate 
practical interest. 
2. Differential equations 
In the first place we have 
© (6") = i n+, (3 
fae} 
o2 ps) 
— (6%) = —&n+? = ——_(€"). (4 
op? O(a?) 


The partial differential equation (4), although simple in appearance, Is 
not very helpful numerically. On the other hand, ordinary differential 
. . . . . . J ] 
equations in 8, which the functions satisfy when « is fixed, are very usefil 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 2 (1957)] 
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asa means of computation. They can be derived from (1) or (2) by integra- 
tion by parts. It can be shown, after some manipulation, that 


Ht 


fe En 2h] 3| n —4 en 4) ip fen 1_£n-3)_ (5) 
| 2 x" | 2a? 2a? 

This linear relation between é-functions of different orders is, in virtue of 
equation (3), an ordinary differential equation in the f-direction of fourth 
order for &"~4, except when n 4, in which case it reduces to a third-order 
equation fo1 6}, 


The real part of equation (5) when n 4 is 
> 


{ d' 9 Cd? -(1 , l \; 5 | n 0 (6) 


; 
\dp3 2° dp* - dp 2x2} 


. 2 
dao 


: being fixed), and numerical integration of this is a very convenient method 
of calculating £.-functions of odd order, and £,-functions of even order. 
The integration may be arranged so as to give a few of the required functions 
simultaneously, others, e.g. HE? and E-', being produced subsequently by 
quadrature of £1 (with due regard to building-up error). Since the second 
derivative of E! (namely, — £%) is available, quadrature formulae more 
rapidly convergent than the usual ones can be used, viz. either 


f = —— ” ~ 
blems 0 | f(x) da hf og + 39°” ise? I (7) 
OTa , 
eans oe lx j he 152 31_54 Te aad "7 
— ol 0 | f(x) da pl Af— 12 pl — 30°44 16800 a | (8) 


in the standard finite-difference notation, for the single quadrature, and 


AQ 
2 { f(a) dade = np 1+ ase—anst..3f" 9 
9 O° | | f(x) dada fe jot 30° 50409 ***f, (9) 
»diat for the double quadrature. 
More accurate values of E>* can, if necessary, be obtained by a separate 
step-by-step integration of, for example, the equation 
B 
- 4 M8 §t—to) . tf os 
E = 2. Eo} 4 E>) dp, (10) 
B B B. 
0 
which follows from (5) with n = 2. In applying (10) the terms in the square 
bracket will be already known. 
Comparable differential equations for the moduli and phases, M,, and 
ce, is >,, exist, but do not seem capable of taking a form nearly as convenient 
ential. for computation as equation (6), or the similar imaginary part of equa- 


tion (5). 
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3. Asymptotic series} 

When f is large enough for asymptotic series to give &” to the required 
accuracy, this becomes a more economic method of computation than 
integration of the differential equations, since values of M,, and ¢,, can 
then be calculated individually and interpolated from a much wider interval 
in 8 than that at which the differential equations can conveniently be 
integrated. In fact M? or M, v8 and ¢,—8 are smooth functions of 1/8 
when f is large, and use of this property allows the whole range of B to be 
covered very concisely. 

One form of asymptotic series, valid for large B/x, can be found by 
expanding the exponential factor in (1). To take the special case n = 1: 

har 


- >) oa — x 
E1+iz! = [ sec {122 see?# + a oe 


48sec 
| 


0 


frag? iM HL (ge_1)-tet8e de 
~ \! x pet ry pit | (: 1) Is. 


1 
x 
Since — | (s?—1)-‘e*°s ds is an integral representation of a Hankel function 
a) 
i 
of zero order, we have 
-_—" 2r 
mI 1. SY * yan 
ta ST “OF ‘a 9 (] 1) 
C 2 Jo r! 0 (B) 
r=0 
— y2" 
El~w 4a S *_Jencp). (12) 
~ hag 
r=0 


Similar series for H-functions of other orders can be found by differentiating 
or integrating (11) and (12) with respect to 8—integrations being naturally 
from infinity. 

Numerically, successive coefficients in these series could be found from 
the recurrence relations which derivatives of Bessel functions satisfy, or 
they could be expressed in terms of Bessel functions of other orders. 

Closely related series can be found in the same way, after removing a 
factor e-**. For example, 

d? 


— - , 0 ee, a 
E!1+4iE! ~ —I4re {ita (1+ J) + 5i(1+ ga +--+} (Yo(B)—th(8)). 
(13) 


are useful 


’ 


+ This term is used to cover series which, while not asymptotic in a strict sense 
for large arguments. 
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In view of the relation 


I 


2 n 
[1+—4) (26a) = (2r—1)(2r—8)...(2r—2n+ Iar-"G,_(@), (14) 
which can be proved by induction, and which holds when @ is either of 
the Bessel functions J or Y, (13) gives 


: ; 1.(8 i amo - 
Ei~w Lo e- YB) 1 q* (P) _, 1. x4 a(P) | (15) 
' | B 2! p? 

J,(B) . 1.3 4d9(B) 
1 = = x } R\_) 2“ 1 5 A Tete 4“2 af: ) 
and E} bar € )70(8) X ~B ay x a i (16) 
Expansions in terms of elementary functions, valid for large B/«?, can 
be obtained at once from these series if the Bessel functions are themselves 

replaced by the well-known asymptotic formulae 


J (B) (2/7B)*{P, (B)cos(B—4naz— }r7)—Q,,(B)sin(8—4nz—}n)} 
Y_(B) (2/7B)*{Q,,(8)cos(B—4n2— }r) +P, (B)sin(B—4na2—}n)} 
in which P(g (4n2 1)(4n? 9) 
b)~ 1— ——____ |. .,, 
2! (88)? 
4n?— ] 


Q),,(B) ~~ - — 


There is a further type of asymptotic expansion of these functions which, 
for numerical purposes, is preferable to the above types, and will therefore 
be treated a little more fully. The first step consists in replacing the 
exponential factor in (1) by an integral representation: 
1 f¢ sia i. 
e—as* e~ “40° cos(su) du. (17) 
XV 7T 


. 


0 


Consider, first, H}; there is no restriction in taking 8 to be positive. Then 


FE} e—**8* sin(Bs)(s?—1)-* ds 

; ; ° 9 1 
e—“" 40° cos(su)sin(sB)(s2?—1)-? duds 

Xx \7T 
i 0 

| n : . 

: ia al sin| (B+-u)s]+-sin| (B—u)s]}(s?—1)-! dsdu, 

0 i 


and since (s?—1)-* sin(as) ds = bardy(x) ifx > 0 (18) 
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and lrJ,(x) ifa< 0, 
x 
EF! — vr [ wider (Bt 4) + J (B—u)t du nie ut/da? J B) d 
14 = 4 € iol P= “o\rP L)5 oa = Jo\U—Pp) au. 
x = a 
0 B Q)+ 
(19)+ 5 
If the quantity {J,(8+u)+4,(8—wu)} is expanded by Taylor’s theoren 
as a power series in w, the series (12) is obtained again. Alternatively, jf 
we use the addition theorem 
= 8 
Jj(B+u) = J(B)Iy(u) +2 ¥ (—)J(B)J,(u), (20 
r=1 
. . . . . oJ ? 
which is valid for all 8 and uw, together with the fact that 
| e-widat J(u) du = avre-“l2] (402), (21 
0 
where /, is the usual modified Bessel function, we find 
; j 
\ -a?/2 { 2 . .: 2\ | \7 2 2 
Es 37 e-% 2) Jo(B)Lo( 40) +2 > J,,(B)L (30°) | i ew 4a" Jy(u—B) du. 
r=1 mX . 
B (22 
The integral in (22) tends rapidly to zero with increasing 8/«, and the series 
alone can conveniently be used to calculate the function for all values of £ 
greater than that at which it first agrees, to the necessary number of 
figures, with the numerical solution of the differential equation. y 
Similar series for the other E,-functions of odd order and £,-functions 
of even order follow by differentiation or integration. Thus, 
C | 2 | 
- 1 aie ares ey , ‘ 
E* — SS E; — 57 e Jo B)Lo(30 ) +2 Zz J9(B)L(4a )j 
Cf r=1 
| ‘ 
— 1p oe-a%2 1,2 12 
= ane I (B)Lo( 30 4 > [Yer-a(B)- Jxy1(B) Go?) |, 
r 
~ M2 mw iee-v2 > 5 2 93 
1.e. E? ~ 47e-% >») So,44(B)bL,.4($0 ), (23 
r=0 
where 6 is the central difference operator in the direction of the order! 
Again ; 


io@) 

3 Ds 12/2 | ; _ 9 | | 
BS ~ Yre-#1Ty—L)— ¥ Soe ps2 + Ip) 
which may be written 

3 
13 1 ve/1 7 S27 4 S27 | ») 
B38 ~ —fre-*2/1, 8+ > dy, d°,), (24 

+ It may be remarked here that numerical quadrature from 0 to o of even functions 

with an ‘e-**’ factor can often be carried out very efficiently by simple summation of the / 


integrand at a surprisingly wide interval, and in (4), where this fact is discussed, the integra 
for E% which corresponds to (19) is treated as an example. 
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where the J-functions are understood to have argument 8 and the J- 
functions to have argument ($a°); 6 is again the central difference operator 
in the r-direction. That this sequence of formulae can be continued is 


perhaps most easily seen by differentiating with respect to a*, using (4) and 


the relation ] 
ui » » 29909 » - 
— fe-@2 7 (4a2)) = fe-*2927 (40°). 25) 
d( : 2 
The sequence continues also in the opposite direction, when central sum- 
mation operators appear in place of central difference operators. It is then 
important from the practical point of view to choose the implied arbitrary 
constants so that the sums tend to zero as r increases. For example we have 
| P in (7 ; | ol 
E | El dB ~ —47e-* *) | J,(8) dB 1($07)+2 5 | 42,(8) dB Lg"); 
. : r 1 2 
which can be reduced to 
= “27 1 l 1 t If 
> ze-% Es I J3(31, T I) T I5(31, 1 I, ] I;) T ..-}}. (=6) 


But the coefficients of J,.,, in (26) do not decrease, and it is better to 


write the series in the form 


7 | | Jo(a) dar) +- 27re-*"S,(B)(3e% ?—31,(3.07)) + Jy(der?—34g—L) 4 


+ J,(fe*21],—I,—I,)+...} (27) 


n which the coefficients 
b1o+ge%?—L—L—...—] (28) 


tend to zero as 7 increases. 
If c denotes the central summation operator in the r-direction (i.e. the 


nverse of 5, previously used), the coefficient (28) may be written —ol,,,. 


iol,., (29) 


r+ 


Thus E® ~ tn 1 ( J, (x) de) Qa e-#2 > J, 


r=0 


U 
where, again, J-functions are understood to have argument f, J-functions 
) 


to have argument («?/2), and 


ol 1 (Iy(4a2)+-e%*2) 


r+3 - U\S 


a a (30) 


By a further integration we have 


E>1= | E°dgB~ | Fedg, 


oc x 
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where F® denotes the right-hand side of (29). That is, 


B : 
E5*~ 4r | (i I(x) de) dt +- 
x 0 
+-2ar € 2h Jy ol, + (Jp +2, )ol; +(Jp+2J,4 2J,)ol;-\-...} 
B r=+s 
{since | J,,.,(2) dx - > J,,(B)) 
snl bin | Jo(x) dx - (8) + 
0 
+ 2a e121 J (ol, +o, 4 ol,+...)4+2,(ol,;+o],+...)+ 
+ 23,(of,+-o1,+...)+...}, 
on reduction of the repeated integral of J. Thus 
E;i~ 18) ( Jy(x) dx 4 (8)| -Arr¢ 2 by 07 ly s Jy, 071}, (31 
0 r=1 
where el, = — s of,.., | 
r=0 : , (32) 
and o* 1. o*I],4 (ol, | oJ, 4 tol, +) | 


which clearly tends to zero as r increases. There are corresponding formulae 
for the £,-functions of odd order and the £,-functions of even order, in 
which Bessel functions of 8 of the second kind, instead of the first, appear. 
The formulae are 


B 
E> ~ 1nB\¥,(B)— | F(x) de t-4re-*2(4¥, o2h+ ¥ Yo,0%,|, (33) 
0 r=1 
B 
E® ~ —ia | Y, (a) da —2me-**2 > Vo,41 01,5, (34) 
0 r=0 
Et ~ —4re-™?/Y,(B)Iy(402) +2 ¥ ¥,(8)1,(302)}, (35) 
| 4 | 
By ~ yre-? Y Yor Shy, (36) 
r=0 
EB? ~ Jre-**{4¥,821,+ > Yo,%F,1. (37) 
r=1 


In these formulae, Y-functions are understood to have argument £, and 
I-functions to have argument (42); § and o have the same meanings as 
before. Although, unlike the J-series, the Y-series do not converge, they 
are very effective for large B/«. 

The accuracy provided by these series depends somewhat on the order 
of the function, as well as upon «, but for the most part four decimals are 
obtainable correctly by about B = 7x. 
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The integrals of Bessel functions occurring in (29), (31), (33), and (34) are 
tabulated at least up to B \0. For larger values of 8 they can be found 
from the asymptotic series 


p 


> (y¥(x) dx = [ ¥,(x) dx 


3 12.3% 1°.3°.5* 1 wunlt 2S, PFS 
~ l — -} B ——————— —_—— 
y, | 32 4 R6 j o(f \B p3 if pe | 
(38) 
| | ( J, (a da [ J,(a ) da 
~ 3. a 
~ —J,(B)i1——+...}+J,(B)} -—_—__—+-...}. (39 
NP) EBT ef TMP Ba To ) 


, 4. Initial values and ascending series 
Expansion of the factor containing f in the original integral leads to 


99) series for 6” in ascending powers of 8 convergent for all values of « and f. 


The convergence is not rapid enough for convenient computation, except 
for small values of B/a, but this is sufficient to enable starting values for 


Lula wil . : ; : : 
integration of the differential equations to be obtained in this way. 
The coefficients in the series, being in fact H-functions with 8 — 0, are 
S 
ea : 
+ all of the form " 
| e~vsec*l seed dé. (40) 
re 0 
They fall naturally into two families according as m is even or odd, corre- 
£ 
sponding to the functions whose asymptotic series contain Bessel functions 
(34) _ of the first and second kinds respectively. Taking, first, m = 2, we have: 
; e-@ 
(35 E?(«, 0) | e—v*sec*? sec? dé Laz ‘ (41) 
P 7 X 
0 
(3¢ , ai 
, From (41) (and (4)) it follows that 
7 dr’ 1 
i 2". 0 r-l1a/y fA-te-A} 42 
E?"(«, 0) yr-t dy qqrai4 a (42) 
ant . 
where A Y and rol 
3 as 
hey \] 0 t -_ e —_ . ’ “ 
Also E(x, 0) | £2 d(a?) = 4n{1—erf(a)}, (43) 
der ¥ 


where 
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and is tabulated. For odd values of m we take first m — 1, and have 
E}(a, 0) = e- sec?) seg 6 dé = te-/2K (402), (44) 
0 


K, being the usual modified Bessel function of the second kind. Then by 
integration and differentiation of (44) with respect to a? there follow: 


Ee (a, 0) = $02 e-%"2{ K, (402) — Ky(40)} (45 
and E2*+l(q.0) = 4(—)r wl fe-Al2K (4 A)}. (46 
c 2 dAr' O\2 j 


A three-term recurrence relation between these coefficients can be obtained 
by setting 8 = Oin (5). The relation is the same for both families, namely, 
n 


E™a, 0) sya EY~3(«, 0) — 4 yn, 0). (47 
| ) y2 


2a? | 2 
Other series of interest arise on further development of formula (22) and 
the formulae derived from (22) by differentiation or integration with 
respect to 8. The difference between EH! and the (convergent) series used 
asymptotically for its calculation is given there as 


_* | e- wide Jy(u—B) du, 


B 
x 
° \i7T ; . B)2/4 2 
1.e. onan a | e+PV ida" J(v) dv (48 
» 
2a . f 
0 
x 
\7 82/42 i 24g? *Bi2a2 
or ee e—P*/4a eV ee eer Ale) ae. (49 
Ae" 





° 8) 
x \ a 
lay : i R\r l si m 
Paine” 340 ee con e-V"da'yr J,(v) dv}, 
2a | 207] r! 
r=0 0 
which may be written 
x 
~ r 
lar Baa? \Y A (a) P x) ; (5( 
/ r a] 
7a) vs 
r=0 
x 
where A,(a) = (—)'*1 | e-Pt'Jy(2at) dt. (51 
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e The first two coefficients in (50) are, in finite terms: 
ro War e- 2] (407), (52) 
4 A, = }e-, (53) 


ind later ones are connected with these by the relations 


= x" on. | (" “|, (54) 
~ d(x?) | 2 


und A, 2 " °)A, (" > +) 4, 2° (55) 


» which can be deduced from the integral (51) by integration by parts. 


\ fairly simple third-order differential equation (in 8 or B/«) exists for 


the function 
me x ‘ 
: ; “ (B/x) — 
(a, 6/a) S A,(«) i —s (56) 
— ,} 
{7 r=0 
but. as a method of calculation, this is inferior to integration of the 
\ani( differential equations given earlier. 
Wl » a 
5. Acknowledgement 
The particular function #%(a,8) has been previously tabulated by the 
Mathematies Division of N.P.L. and thanks are due to Dr. E. T. Goodwin 
ind Mr. F. W. J. Olver of that Division for making available their investi- 
gations on the subject. In particular, formulae of the same type as (11), (19). 
‘6), and (47), and the asymptotic expansion in terms of elementary 
functions are due to them. 
ts 
REFERENCES 
1. T. H. Havetock, Proc. Roy. Soc. A, 108 (1925) 582. 
2. Quart. J. Mech. App. Math. 5 (1952) 129. 
{ 3. H. L. Ponn,. David Taylor Model Basin, teport 795 (Nov. 1951). 
4. E. T. Goopwry, Proc. Camb. Phil. Soc. 45 (1949) 241. 

















ON THE GENERALIZED BER, BEI, KER, AND 
KET FUNCTIONS WITH APPLICATION TO PLATE 
PROBLEMS 


By YI-YUAN YU (Syracuse University, Syracuse) 


[Received 24 May 1956] 


SUMMARY 

It is pointed out in this note that the solution of the differential equation 
V2V2y—sV2y+ty 0 when s?—4t < 0 
may best be put in terms of the real and imaginary parts of the modified Bess: 
functions of the zero.,order, denoted by J)(ax),, [,(ax);, Ko(ax),, and K,(ax);. Sine 
these reduce to the ber, bei, ker, and kei functions when s is zero, they are called 
the generalized ber, bei, ker, and kei functions. Their application to plate problems 
is indicated. 


1. Generalized ber, bei, ker, kei functions 
IN certain plate problems it is necessary to solve the differential equation 
V2V2y—sV2y+ty = 0, (1 
J eo ae 
where si -- ; 
Ga" 2éz 
x, y are real variables, and s, ¢ are real constants. When s?— 4¢ is positive 
or zero, the solution of (1) may be conveniently put in terms of the ordinary 
Bessel functions. For instance, when both s?—4t and s itself are positive, 
the solution is of the form 


y = CO, [y(ax)+C, Iy(Bax)+-C; Ko(ax)+C, Ko(Bx), (2) 
s+ (s?—4t)*]4 
where x, B= | _ste =r lt (3) 


However, when s?— 4 is negative, « and 8B become complex conjugates to 
each other, and so are J,(«a) and J,(8x), and K,(aa) and K,(f2), all of which 
now have complex arguments. It is then more convenient to write (2) in 
the form 
y = C, Lax), +C, [y(aa);+C; Ko(ax),+-C, Ko(ar);, (4) 
where the subscripts 7 and 7 indicate the pure real and pure imaginary parts, 
respectively, of the function to which they are attached. 
The four functions in (4) reduce to the ber, bei, ker, and kei functions 
when s is zero. Again, we see that the same is true when the real and 
imaginary parts of the complex number « are equal to each other. These 
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functions are therefore called the generalized ber, bei, ker, kei functions. 
r) They have similar continuity properties as the ber, bei, ker, and kei func- 
tions, namely, J,(z), and J,(z); are continuous for all real values of z, and so 
are K,(z), and K,(z); except when z 0. 
, Similar functions of an order v other than zero may also be formulated, 
such as [,(z),, [(z);, ete. 


The following formulae for the derivatives are easily obtained: 


d 
I, (x2), al,( vx), —bI, ( XX) 55 


dx 
1 Lp \ xX); al,(ax);+bL,(ax),, 
( l 
essel | j I vx), al,( YX). bl, ( XX) ; L( ut)... 
o. aa x 


" l ] 
ble L YX); al, (aa); + b1,(a2x),- [,(a2);, 


da x 
—_ ; : 
-Ko(ax), ak,(aa),+bK,(ax),, 
Ato! d ‘ » , 
j K,( XX); ak,( XX); bK,( XX) ps 
; : Ds 
; K, (a2), ak,(ax),+bK,(ax);- K,(ax),, 
ax x 
itive d , » , ] y 
Ky (ax); ak, (aa), —bK,(ax),— — Ky (ax);, 
nary AX x 


tive, where the constants a, b, a,, and b, are defined by 
a+b Y, a,+-0b, x”, 


> These formulae are useful in solving certain thin plate problems. 
2. Application to plate problems 


We shall consider the problem of axially symmetrical bending or buckling 
ofa thin plate (flexural rigidity D) on elastic foundation (modulus k) under 


+ ) the simultaneous action of lateral pressure p and force in the middle plane NV 
“7 2 (tension positive). According to the classical theory, the deflexion of the 
. plate is governed by the equation 
DV?V*y—N V2y+ky = p. 
— Since the left-hand member of this equation is of the same form as that 
of (1), its complementary function is also given by (4), in which the value 
ons} oof wis now 


N s 


and N- i(4kD- N?)*]4 
hese 2D ; 
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Corresponding to this complementary function for the deflexion y, the 
slope, radial bending moment, and shearing force may be calculated straight. 
forwardly by means of the differentiation formulae given above; thus, 
dy 


== = C,[al,(ax),—bL, (ax); |+C,[al, (ax); +bL,(ax),|— 
dx 


—C,[ aK, (a2),—bK, (ax), |—C,[aK, (ax); +-bK,(aa),], 


l?7y | vdy 
|| —s ant “a i 
I, , (4 x 2) 


= — plc fer Ja vv),—b, Ih(ax); — : lah vr),—bI, (ax 


C4; I,(ax);+b, Iy(axx),— ; ~ ” [al,( var); + bE, (a2), 


C,1a, Ko(ox),—b, Ka(ax), +=" [aK,(ax),—bK,( 
x 


Cs) r 
; it, 
7 Cy Ko( xX); +b, Ko( YX), x r : [aky( wx) 5 


d _dy 
) V2y)+N — 
dx \ y) dx 


D[Ci{(a, a+b, b) 1, (aa),—(a, b—b, a) I, (ax);} 
-O,f{(a, b—b, a) I, (ax),+ (a, a+b, b)I,(ax);}— 
—C,{(a, a+b, b) Ky (ax),— (a, b—b, a) Ky (ax);}— 


a) 
C,{(a, b—b, a) Ky (ax),+ (a, a+b, b)K,(ax);}]. 
Any axially symmetrical bending or buckling problem of a thin plate may 
then be solved in the usual manner. 
On the basis of Reissner’s plate theory, Naghdi and Rowley (1) derived 
the following equation for the axially symmetrical bending of plates om 
elastic foundation, without any force in the middle plane: 


V2V2L—BVL+C = 0, 


which is obviously a special case of (1). The solution may therefore be put 
in the same closed form as (4). These authors, however, used the ber, bei, 
ker, kei functions together with series in powers of 6?, which is rather 
complicated. 
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